/* ------------------------------------------------------------- file: versteeg.cpp This is a simple test of PETSC Matrix operations and Linear equation solvers. This problem comes from Example 7.2 in Versteeg, H.K. and W. Malalasekera, 1995. An introduction to computational fluid dynamics, the finite volume method. Prentice Hall. 257 pp. ------------------------------------------------------------- */ /* ------------------------------------------------------------- Created October 10, 2002 by William A. Perkins Last Change: 2013-09-26 14:32:33 d3g096 ------------------------------------------------------------- */ static const char* SCCS_ID = "$Id: versteeg.c,v 1.3 2008/05/24 23:04:39 d3g096 Exp $ Battelle PNNL"; static char help[] = "Solves an example heat transfer problem from Versteeg and Malalasekera with SLES.\n\ The following options are available:\n\ -view_solution yes/no: write the initial and final solution to stdout\n\ -problem yes/no: cause problem by choosing local ownership\n\n"; #include #include /* ------------------------------------------------------------- assemble ------------------------------------------------------------- */ int assemble(const int imax, const int jmax, const double scale, Mat A, Vec b) { static const float k = 1000; /* conductivity, W/m/K */ static const float t = 0.01; /* plate thickness, m */ static const float W = 0.3; /* plate width, m */ static const float H = 0.4; /* plate height, m */ PetscErrorCode ierr; int i, j; float ap, aw, ae, as, an, bp; float dx = W/(float)imax; float dy = H/(float)jmax; int iP, iN, iS, iE, iW; PetscScalar v; PetscInt lo, hi; ierr = MatGetOwnershipRange(A, &lo, &hi); CHKERRQ(ierr); for (i = 0; i < imax; i++) { for (j = 0; j < jmax; j++) { iP = i*jmax + j; if (! (lo <= iP && iP < hi) ) continue; iE = (i+1)*jmax + j; iW = (i-1)*jmax + j; iN = i*jmax + (j+1); iS = i*jmax + (j-1); bp = 0.0; ap = 0.0; if (j == 0) { /* insulated south boundary */ as = 0.0; bp += 0.0; ap -= 0.0; } else { as = (k/dx)*(dx*t); } if (j == jmax - 1) { /* constant tempurature (100C) north boundary */ an = 0.0; bp += 2*k/dy*(dy*t)*100.0; ap -= -2*k/dy*(dy*t); } else { an = (k/dx)*(dx*t); } if (i == 0) { /* constant flux (500kw/m2) west boundary */ aw = 0.0; bp += 500000.0*dy*t; ap -= 0.0; } else { aw = (k/dx)*(dx*t); } if (i == imax - 1) { /* insulated east boundary */ ae = 0.0; bp += 0.0; ap -= 0.0; } else { ae = (k/dx)*(dx*t); } ap += as + an + aw + ae; v = ap*scale; ierr = MatSetValues(A,1,&iP,1,&iP,&v,INSERT_VALUES); CHKERRQ(ierr); if (an != 0.0) { v = -an*scale; ierr = MatSetValues(A,1,&iP,1,&iN,&v,INSERT_VALUES); CHKERRQ(ierr); } if (as != 0.0) { v = -as*scale; ierr = MatSetValues(A,1,&iP,1,&iS,&v,INSERT_VALUES); CHKERRQ(ierr); } if (ae != 0.0) { v = -ae*scale; ierr = MatSetValues(A,1,&iP,1,&iE,&v,INSERT_VALUES); CHKERRQ(ierr); } if (aw != 0.0) { v = -aw*scale; ierr = MatSetValues(A,1,&iP,1,&iW,&v,INSERT_VALUES); CHKERRQ(ierr); } v = bp*scale; ierr = VecSetValues(b, 1, &iP, &v, INSERT_VALUES); CHKERRQ(ierr); } } ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); return(0); } /* ------------------------------------------------------------- Main Program ------------------------------------------------------------- */ int main(int argc, char **args) { static const int imax = 3; static const int jmax = 4; PetscBool print_solution(PETSC_FALSE); PetscBool set_local_size(PETSC_FALSE); int its; PetscErrorCode ierr; Vec x, b; Mat A; PetscScalar v; KSP ksp; VecScatter ctx; int rank; int nproc; PetscInt local_size(PETSC_DECIDE); PetscInt global_size(imax*jmax); PetscInitialize(&argc,&args,(char *)0,help); PetscOptionsGetBool(NULL, "-print_solution", &print_solution, NULL); PetscOptionsGetBool(NULL, "-problem", &set_local_size, NULL); ierr = MPI_Comm_size(PETSC_COMM_WORLD,&nproc); ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); ierr = PetscSplitOwnership(PETSC_COMM_WORLD, &local_size, &global_size); CHKERRQ(ierr); if (set_local_size) { // If there is more than one processor, alter the default local // ownership so that process 0 has one less row and the last // process has more row. if (nproc > 1) { if (rank == 0) { local_size -= 1; } else if (rank == nproc - 1) { local_size += 1; } } } ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRQ(ierr); ierr = MatSetSizes(A, local_size, PETSC_DECIDE, imax*jmax, imax*jmax); CHKERRXX(ierr); ierr = MatSetFromOptions(A); CHKERRQ(ierr); ierr = MatSetUp(A); CHKERRQ(ierr); ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); ierr = VecSetSizes(x,local_size,imax*jmax);CHKERRQ(ierr); ierr = VecSetFromOptions(x);CHKERRQ(ierr); ierr = VecDuplicate(x,&b);CHKERRQ(ierr); ierr = assemble(imax, jmax, 1.0, A, b); CHKERRQ(ierr); v = 0.0; ierr = VecSet(x,v);CHKERRQ(ierr); if (print_solution) { ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); } ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE); CHKERRQ(ierr); ierr = KSPSetTolerances(ksp, 1e-06, 1e-12, PETSC_DEFAULT, PETSC_DEFAULT); CHKERRQ(ierr); ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); ierr = KSPSetUp(ksp); CHKERRQ(ierr); ierr = KSPSolve(ksp, b, x);CHKERRQ(ierr); ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"iterations %d\n",its);CHKERRQ(ierr); if (print_solution) { ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); } ierr = KSPDestroy(&ksp);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = PetscFinalize(); return(ierr); }