Actual source code: ex35.cxx

      1: /*
      2: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation

      4:    -div \rho grad u = f,  0 < x,y < 1,

      6: Problem 1: (Default)

      8:   Use the following exact solution with Dirichlet boundary condition

     10:     u = sin(pi*x)*sin(pi*y)

     12:   and generate an appropriate forcing function to measure convergence.

     14: Usage:

     16:   Measure convergence rate with uniform refinement with the options: "-problem 1 -error".

     18:     mpiexec -n $NP ./ex35 -problem 1 -error -n 16 -levels 5 -pc_type gamg
     19:     mpiexec -n $NP ./ex35 -problem 1 -error -n 32 -levels 4 -pc_type gamg
     20:     mpiexec -n $NP ./ex35 -problem 1 -error -n 64 -levels 3 -pc_type mg
     21:     mpiexec -n $NP ./ex35 -problem 1 -error -n 128 -levels 2 -pc_type mg
     22:     mpiexec -n $NP ./ex35 -problem 1 -error -n 256 -levels 1 -mg
     23:     mpiexec -n $NP ./ex35 -problem 1 -error -n 512 -levels 0 -mg

     25:   Or with an external mesh file representing [0, 1]^2,

     27:     mpiexec -n $NP ./ex35 -problem 1 -file ./external_mesh.h5m -levels 2 -error -pc_type mg

     29: Problem 2:

     31:   Use the forcing function

     33:      f = e^{-((x-xr)^2+(y-yr)^2)/\nu}

     35:   with Dirichlet boundary conditions

     37:      u = f(x,y) for x = 0, x = 1, y = 0, y = 1

     39:   or pure Neumman boundary conditions

     41: Usage:

     43:   Run with different values of \rho and \nu (problem 1) to control diffusion and gaussian source spread. This uses the internal mesh generator implemented in DMMoab.

     45:     mpiexec -n $NP ./ex35 -problem 2 -n 20 -nu 0.02 -rho 0.01
     46:     mpiexec -n $NP ./ex35 -problem 2 -n 40 -nu 0.01 -rho 0.005 -io -ksp_monitor
     47:     mpiexec -n $NP ./ex35 -problem 2 -n 80 -nu 0.01 -rho 0.005 -io -ksp_monitor -pc_type hypre
     48:     mpiexec -n $NP ./ex35 -problem 2 -n 160 -bc neumann -nu 0.005 -rho 0.01 -io
     49:     mpiexec -n $NP ./ex35 -problem 2 -n 320 -bc neumann -nu 0.001 -rho 1 -io

     51:   Or with an external mesh file representing [0, 1]^2,

     53:     mpiexec -n $NP ./ex35 -problem 2 -file ./external_mesh.h5m -levels 1 -pc_type gamg

     55: Problem 3:

     57:   Use the forcing function

     59:      f = \nu*sin(\pi*x/LX)*sin(\pi*y/LY)

     61:   with Dirichlet boundary conditions

     63:      u = 0.0, for x = 0, x = 1, y = 0, y = 1 (outer boundary) &&
     64:      u = 1.0, for (x-0.5)^2 + (y-0.5)^2 = 0.01

     66: Usage:

     68:   Now, you could alternately load an external MOAB mesh that discretizes the unit square and use that to run the solver.

     70:     mpiexec -n $NP ./ex35 -problem 3 -file input/square_with_hole.h5m -mg
     71:     mpiexec -n $NP ./ex35 -problem 3 -file input/square_with_hole.h5m -mg -levels 2 -io -ksp_monitor
     72:     mpiexec -n $NP ./ex35 -problem 3 -file input/square_with_hole.h5m -io -ksp_monitor -pc_type hypre

     74: */

     76: static char help[] = "\
     77:                       Solves the 2D inhomogeneous Laplacian equation.\n \
     78:                       Usage: ./ex35 -problem 1 -error -n 4 -levels 3 -mg\n  \
     79:                       Usage: ./ex35 -problem 2 -n 80 -nu 0.01 -rho 0.005 -io -ksp_monitor -pc_type gamg\n  \
     80:                       Usage: ./ex35 -problem 3 -file input/square_with_hole.h5m -mg\n";

     82: /* PETSc includes */
     83: #include <petscksp.h>
     84: #include <petscdmmoab.h>

     86: #define LOCAL_ASSEMBLY

     88: typedef enum {
     89:   DIRICHLET,
     90:   NEUMANN
     91: } BCType;

     93: typedef struct {
     94:   PetscInt  dim, n, problem, nlevels;
     95:   PetscReal rho;
     96:   PetscReal bounds[6];
     97:   PetscReal xref, yref;
     98:   PetscReal nu;
     99:   PetscInt  VPERE;
    100:   BCType    bcType;
    101:   PetscBool use_extfile, io, error, usetri, usemg;
    102:   char      filename[PETSC_MAX_PATH_LEN];
    103: } UserContext;

    105: static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
    106: static PetscErrorCode ComputeRHS(KSP, Vec, void *);
    107: static PetscErrorCode ComputeDiscreteL2Error(KSP, Vec, UserContext *);
    108: static PetscErrorCode InitializeOptions(UserContext *);

    110: int main(int argc, char **argv)
    111: {
    112:   KSP ksp;
    113:   PC  pc;
    114:   Mat R;
    115:   DM  dm, dmref, *dmhierarchy;

    117:   UserContext user;
    118:   const char *fields[1] = {"T-Variable"};
    119:   PetscInt    k;
    120:   Vec         b, x, errv;

    122:   PetscFunctionBeginUser;
    123:   PetscCall(PetscInitialize(&argc, &argv, nullptr, help));

    125:   PetscCall(InitializeOptions(&user));

    127:   /* Create the DM object from either a mesh file or from in-memory structured grid */
    128:   if (user.use_extfile) {
    129:     PetscCall(DMMoabLoadFromFile(PETSC_COMM_WORLD, user.dim, 1, user.filename, "", &dm));
    130:   } else {
    131:     PetscCall(DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, user.usetri, user.bounds, user.n, 1, &dm));
    132:   }
    133:   PetscCall(DMSetFromOptions(dm));
    134:   PetscCall(DMMoabSetFieldNames(dm, 1, fields));

    136:   /* SetUp the data structures for DMMOAB */
    137:   PetscCall(DMSetUp(dm));

    139:   PetscCall(DMSetApplicationContext(dm, &user));

    141:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    142:   PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, &user));
    143:   PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, &user));

    145:   if (user.nlevels) {
    146:     PetscCall(KSPGetPC(ksp, &pc));
    147:     PetscCall(PetscMalloc(sizeof(DM) * (user.nlevels + 1), &dmhierarchy));
    148:     for (k = 0; k <= user.nlevels; k++) dmhierarchy[k] = NULL;

    150:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of mesh hierarchy levels: %d\n", user.nlevels));
    151:     PetscCall(DMMoabGenerateHierarchy(dm, user.nlevels, NULL));

    153:     /* coarsest grid = 0, finest grid = nlevels */
    154:     dmhierarchy[0]         = dm;
    155:     PetscBool usehierarchy = PETSC_FALSE;
    156:     if (usehierarchy) {
    157:       PetscCall(DMRefineHierarchy(dm, user.nlevels, &dmhierarchy[1]));
    158:     } else {
    159:       for (k = 1; k <= user.nlevels; k++) PetscCall(DMRefine(dmhierarchy[k - 1], MPI_COMM_NULL, &dmhierarchy[k]));
    160:     }
    161:     dmref = dmhierarchy[user.nlevels];
    162:     PetscCall(PetscObjectReference((PetscObject)dmref));

    164:     if (user.usemg) {
    165:       PetscCall(PCSetType(pc, PCMG));
    166:       PetscCall(PCMGSetLevels(pc, user.nlevels + 1, NULL));
    167:       PetscCall(PCMGSetType(pc, PC_MG_MULTIPLICATIVE));
    168:       PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
    169:       PetscCall(PCMGSetCycleType(pc, PC_MG_CYCLE_V));
    170:       PetscCall(PCMGSetNumberSmooth(pc, 2));

    172:       for (k = 1; k <= user.nlevels; k++) {
    173:         PetscCall(DMCreateInterpolation(dmhierarchy[k - 1], dmhierarchy[k], &R, NULL));
    174:         PetscCall(PCMGSetInterpolation(pc, k, R));
    175:         PetscCall(MatDestroy(&R));
    176:       }
    177:     }

    179:     for (k = 1; k <= user.nlevels; k++) PetscCall(DMDestroy(&dmhierarchy[k]));
    180:     PetscCall(PetscFree(dmhierarchy));
    181:   } else {
    182:     dmref = dm;
    183:     PetscCall(PetscObjectReference((PetscObject)dm));
    184:   }

    186:   PetscCall(KSPSetDM(ksp, dmref));
    187:   PetscCall(KSPSetFromOptions(ksp));

    189:   /* Perform the actual solve */
    190:   PetscCall(KSPSolve(ksp, NULL, NULL));
    191:   PetscCall(KSPGetSolution(ksp, &x));
    192:   PetscCall(KSPGetRhs(ksp, &b));

    194:   if (user.error) {
    195:     PetscCall(VecDuplicate(b, &errv));
    196:     PetscCall(ComputeDiscreteL2Error(ksp, errv, &user));
    197:     PetscCall(VecDestroy(&errv));
    198:   }

    200:   if (user.io) {
    201:     /* Write out the solution along with the mesh */
    202:     PetscCall(DMMoabSetGlobalFieldVector(dmref, x));
    203: #ifdef MOAB_HAVE_HDF5
    204:     PetscCall(DMMoabOutput(dmref, "ex35.h5m", NULL));
    205: #else
    206:     /* MOAB does not support true parallel writers that aren't HDF5 based
    207:        And so if you are using VTK as the output format in parallel,
    208:        the data could be jumbled due to the order in which the processors
    209:        write out their parts of the mesh and solution tags */
    210:     PetscCall(DMMoabOutput(dmref, "ex35.vtk", NULL));
    211: #endif
    212:   }

    214:   /* Cleanup objects */
    215:   PetscCall(KSPDestroy(&ksp));
    216:   PetscCall(DMDestroy(&dmref));
    217:   PetscCall(DMDestroy(&dm));
    218:   PetscCall(PetscFinalize());
    219:   return 0;
    220: }

    222: PetscReal ComputeDiffusionCoefficient(PetscReal coords[3], UserContext *user)
    223: {
    224:   switch (user->problem) {
    225:   case 2:
    226:     if ((coords[0] > user->bounds[1] / 3.0) && (coords[0] < 2.0 * user->bounds[1] / 3.0) && (coords[1] > user->bounds[3] / 3.0) && (coords[1] < 2.0 * user->bounds[3] / 3.0)) {
    227:       return user->rho;
    228:     } else {
    229:       return 1.0;
    230:     }
    231:   case 1:
    232:   case 3:
    233:   default:
    234:     return user->rho;
    235:   }
    236: }

    238: PetscScalar ExactSolution(PetscReal coords[3], UserContext *user)
    239: {
    240:   switch (user->problem) {
    241:   case 1:
    242:     return sin(PETSC_PI * coords[0] / user->bounds[1]) * sin(PETSC_PI * coords[1] / user->bounds[3]);
    243:   case 3:
    244:   case 2:
    245:   default:
    246:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Exact solution for -problem = [%" PetscInt_FMT "] is not available.", user->problem);
    247:   }
    248: }

    250: PetscScalar ComputeForcingFunction(PetscReal coords[3], UserContext *user)
    251: {
    252:   switch (user->problem) {
    253:   case 3:
    254:     return user->nu * sin(PETSC_PI * coords[0] / user->bounds[1]) * sin(PETSC_PI * coords[1] / user->bounds[3]);
    255:   case 2:
    256:     return PetscExpScalar(-((coords[0] - user->xref) * (coords[0] - user->xref) + (coords[1] - user->yref) * (coords[1] - user->yref)) / user->nu);
    257:   case 1:
    258:   default:
    259:     return PETSC_PI * PETSC_PI * ComputeDiffusionCoefficient(coords, user) * (1.0 / user->bounds[1] / user->bounds[1] + 1.0 / user->bounds[3] / user->bounds[3]) * sin(PETSC_PI * coords[0] / user->bounds[1]) * sin(PETSC_PI * coords[1] / user->bounds[3]);
    260:   }
    261: }

    263: #define BCHECKEPS                   1e-10
    264: #define BCHECK(coordxyz, truetrace) (coordxyz < truetrace + BCHECKEPS && coordxyz > truetrace - BCHECKEPS)

    266: PetscScalar EvaluateStrongDirichletCondition(PetscReal coords[3], UserContext *user)
    267: {
    268:   switch (user->problem) {
    269:   case 3:
    270:     if (BCHECK(coords[0], user->bounds[0]) || BCHECK(coords[0], user->bounds[1]) || BCHECK(coords[1], user->bounds[2]) || BCHECK(coords[1], user->bounds[3])) return 0.0;
    271:     else // ( coords[0]*coords[0] + coords[1]*coords[1] < 0.04 + BCHECKEPS)
    272:       return 1.0;
    273:   case 2:
    274:     return ComputeForcingFunction(coords, user);
    275:   case 1:
    276:   default:
    277:     return ExactSolution(coords, user);
    278:   }
    279: }

    281: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ptr)
    282: {
    283:   UserContext              *user = (UserContext *)ptr;
    284:   DM                        dm;
    285:   PetscInt                  dof_indices[4];
    286:   PetscBool                 dbdry[4];
    287:   PetscReal                 vpos[4 * 3];
    288:   PetscScalar               ff;
    289:   PetscInt                  i, q, nconn, nc, npoints;
    290:   const moab::EntityHandle *connect;
    291:   const moab::Range        *elocal;
    292:   moab::Interface          *mbImpl;
    293:   PetscScalar               localv[4];
    294:   PetscReal                *phi, *phypts, *jxw;
    295:   PetscBool                 elem_on_boundary;
    296:   PetscQuadrature           quadratureObj;

    298:   PetscFunctionBegin;
    299:   PetscCall(KSPGetDM(ksp, &dm));

    301:   /* reset the RHS */
    302:   PetscCall(VecSet(b, 0.0));

    304:   PetscCall(DMMoabFEMCreateQuadratureDefault(2, user->VPERE, &quadratureObj));
    305:   PetscCall(PetscQuadratureGetData(quadratureObj, NULL, &nc, &npoints, NULL, NULL));
    306:   PetscCall(PetscMalloc3(user->VPERE * npoints, &phi, npoints * 3, &phypts, npoints, &jxw));

    308:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
    309:   PetscCall(DMMoabGetInterface(dm, &mbImpl));
    310:   PetscCall(DMMoabGetLocalElements(dm, &elocal));

    312:   /* loop over local elements */
    313:   for (moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
    314:     const moab::EntityHandle ehandle = *iter;

    316:     /* Get connectivity information: */
    317:     PetscCall(DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect));
    318:     PetscCheck(nconn == 3 || nconn == 4, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only TRI3/QUAD4 element bases are supported in the current example. n(Connectivity)=%" PetscInt_FMT ".", nconn);

    320:     PetscCall(PetscArrayzero(localv, nconn));

    322:     /* get the coordinates of the element vertices */
    323:     PetscCall(DMMoabGetVertexCoordinates(dm, nconn, connect, vpos));

    325:     /* get the local DoF numbers to appropriately set the element contribution in the operator */
    326: #ifdef LOCAL_ASSEMBLY
    327:     PetscCall(DMMoabGetFieldDofsLocal(dm, nconn, connect, 0, dof_indices));
    328: #else
    329:     PetscCall(DMMoabGetFieldDofs(dm, nconn, connect, 0, dof_indices));
    330: #endif

    332:     /* 1) compute the basis functions and the derivatives wrt x and y directions
    333:        2) compute the quadrature points transformed to the physical space */
    334:     PetscCall(DMMoabFEMComputeBasis(2, nconn, vpos, quadratureObj, phypts, jxw, phi, NULL));

    336:     /* Compute function over the locally owned part of the grid */
    337:     for (q = 0; q < npoints; ++q) {
    338:       ff = ComputeForcingFunction(&phypts[3 * q], user);

    340:       for (i = 0; i < nconn; ++i) localv[i] += jxw[q] * phi[q * nconn + i] * ff;
    341:     }

    343:     /* check if element is on the boundary */
    344:     PetscCall(DMMoabIsEntityOnBoundary(dm, ehandle, &elem_on_boundary));

    346:     /* apply dirichlet boundary conditions */
    347:     if (elem_on_boundary && user->bcType == DIRICHLET) {
    348:       /* get the list of nodes on boundary so that we can enforce dirichlet conditions strongly */
    349:       PetscCall(DMMoabCheckBoundaryVertices(dm, nconn, connect, dbdry));

    351:       for (i = 0; i < nconn; ++i) {
    352:         if (dbdry[i]) { /* dirichlet node */
    353:           /* think about strongly imposing dirichlet */
    354:           localv[i] = EvaluateStrongDirichletCondition(&vpos[3 * i], user);
    355:         }
    356:       }
    357:     }

    359: #ifdef LOCAL_ASSEMBLY
    360:     /* set the values directly into appropriate locations. Can alternately use VecSetValues */
    361:     PetscCall(VecSetValuesLocal(b, nconn, dof_indices, localv, ADD_VALUES));
    362: #else
    363:     PetscCall(VecSetValues(b, nconn, dof_indices, localv, ADD_VALUES));
    364: #endif
    365:   }

    367:   /* force right-hand side to be consistent for singular matrix */
    368:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
    369:   if (user->bcType == NEUMANN) {
    370:     MatNullSpace nullspace;
    371:     PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
    372:     PetscCall(MatNullSpaceRemove(nullspace, b));
    373:     PetscCall(MatNullSpaceDestroy(&nullspace));
    374:   }

    376:   /* Restore vectors */
    377:   PetscCall(VecAssemblyBegin(b));
    378:   PetscCall(VecAssemblyEnd(b));
    379:   PetscCall(PetscFree3(phi, phypts, jxw));
    380:   PetscCall(PetscQuadratureDestroy(&quadratureObj));
    381:   PetscFunctionReturn(PETSC_SUCCESS);
    382: }

    384: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, PetscCtx ctx)
    385: {
    386:   UserContext              *user = (UserContext *)ctx;
    387:   DM                        dm;
    388:   PetscInt                  i, j, q, nconn, nglobale, nglobalv, nc, npoints, hlevel;
    389:   PetscInt                  dof_indices[4];
    390:   PetscReal                 vpos[4 * 3], rho;
    391:   PetscBool                 dbdry[4];
    392:   const moab::EntityHandle *connect;
    393:   const moab::Range        *elocal;
    394:   moab::Interface          *mbImpl;
    395:   PetscBool                 elem_on_boundary;
    396:   PetscScalar               array[4 * 4];
    397:   PetscReal                *phi, *dphi[2], *phypts, *jxw;
    398:   PetscQuadrature           quadratureObj;

    400:   PetscFunctionBeginUser;
    401:   PetscCall(KSPGetDM(ksp, &dm));

    403:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
    404:   PetscCall(DMMoabGetInterface(dm, &mbImpl));
    405:   PetscCall(DMMoabGetLocalElements(dm, &elocal));
    406:   PetscCall(DMMoabGetSize(dm, &nglobale, &nglobalv));
    407:   PetscCall(DMMoabGetHierarchyLevel(dm, &hlevel));
    408:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ComputeMatrix: Level = %d, N(elements) = %d, N(vertices) = %d \n", hlevel, nglobale, nglobalv));

    410:   PetscCall(DMMoabFEMCreateQuadratureDefault(2, user->VPERE, &quadratureObj));
    411:   PetscCall(PetscQuadratureGetData(quadratureObj, NULL, &nc, &npoints, NULL, NULL));
    412:   PetscCall(PetscMalloc5(user->VPERE * npoints, &phi, user->VPERE * npoints, &dphi[0], user->VPERE * npoints, &dphi[1], npoints * 3, &phypts, npoints, &jxw));

    414:   /* loop over local elements */
    415:   for (moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
    416:     const moab::EntityHandle ehandle = *iter;

    418:     // Get connectivity information:
    419:     PetscCall(DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect));
    420:     PetscCheck(nconn == 3 || nconn == 4, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only QUAD4 or TRI3 element bases are supported in the current example. Connectivity=%" PetscInt_FMT ".", nconn);

    422:     /* compute the mid-point of the element and use a 1-point lumped quadrature */
    423:     PetscCall(DMMoabGetVertexCoordinates(dm, nconn, connect, vpos));

    425:     /* get the global DOF number to appropriately set the element contribution in the RHS vector */
    426: #ifdef LOCAL_ASSEMBLY
    427:     PetscCall(DMMoabGetFieldDofsLocal(dm, nconn, connect, 0, dof_indices));
    428: #else
    429:     PetscCall(DMMoabGetFieldDofs(dm, nconn, connect, 0, dof_indices));
    430: #endif

    432:     /* 1) compute the basis functions and the derivatives wrt x and y directions
    433:        2) compute the quadrature points transformed to the physical space */
    434:     PetscCall(DMMoabFEMComputeBasis(2, nconn, vpos, quadratureObj, phypts, jxw, phi, dphi));

    436:     PetscCall(PetscArrayzero(array, nconn * nconn));

    438:     /* Compute function over the locally owned part of the grid */
    439:     for (q = 0; q < npoints; ++q) {
    440:       /* compute the inhomogeneous (piece-wise constant) diffusion coefficient at the quadrature point
    441:         -- for large spatial variations (within an element), embed this property evaluation inside the quadrature loop
    442:       */
    443:       rho = ComputeDiffusionCoefficient(&phypts[q * 3], user);

    445:       for (i = 0; i < nconn; ++i) {
    446:         for (j = 0; j < nconn; ++j) array[i * nconn + j] += jxw[q] * rho * (dphi[0][q * nconn + i] * dphi[0][q * nconn + j] + dphi[1][q * nconn + i] * dphi[1][q * nconn + j]);
    447:       }
    448:     }

    450:     /* check if element is on the boundary */
    451:     PetscCall(DMMoabIsEntityOnBoundary(dm, ehandle, &elem_on_boundary));

    453:     /* apply dirichlet boundary conditions */
    454:     if (elem_on_boundary && user->bcType == DIRICHLET) {
    455:       /* get the list of nodes on boundary so that we can enforce dirichlet conditions strongly */
    456:       PetscCall(DMMoabCheckBoundaryVertices(dm, nconn, connect, dbdry));

    458:       for (i = 0; i < nconn; ++i) {
    459:         if (dbdry[i]) { /* dirichlet node */
    460:           /* think about strongly imposing dirichlet */
    461:           for (j = 0; j < nconn; ++j) {
    462:             /* TODO: symmetrize the system - need the RHS */
    463:             array[i * nconn + j] = 0.0;
    464:           }
    465:           array[i * nconn + i] = 1.0;
    466:         }
    467:       }
    468:     }

    470:     /* set the values directly into appropriate locations. */
    471: #ifdef LOCAL_ASSEMBLY
    472:     PetscCall(MatSetValuesLocal(jac, nconn, dof_indices, nconn, dof_indices, array, ADD_VALUES));
    473: #else
    474:     PetscCall(MatSetValues(jac, nconn, dof_indices, nconn, dof_indices, array, ADD_VALUES));
    475: #endif
    476:   }

    478:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    479:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));

    481:   if (user->bcType == NEUMANN) {
    482:     MatNullSpace nullspace;

    484:     PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
    485:     PetscCall(MatSetNullSpace(J, nullspace));
    486:     PetscCall(MatNullSpaceDestroy(&nullspace));
    487:   }
    488:   PetscCall(PetscFree5(phi, dphi[0], dphi[1], phypts, jxw));
    489:   PetscCall(PetscQuadratureDestroy(&quadratureObj));
    490:   PetscFunctionReturn(PETSC_SUCCESS);
    491: }

    493: PetscErrorCode ComputeDiscreteL2Error(KSP ksp, Vec err, UserContext *user)
    494: {
    495:   DM                 dm;
    496:   Vec                sol;
    497:   PetscReal          vpos[3];
    498:   const PetscScalar *x;
    499:   PetscScalar       *e;
    500:   PetscReal          l2err = 0.0, linferr = 0.0, global_l2, global_linf;
    501:   PetscInt           dof_index, N;
    502:   const moab::Range *ownedvtx;

    504:   PetscFunctionBegin;
    505:   PetscCall(KSPGetDM(ksp, &dm));

    507:   /* get the solution vector */
    508:   PetscCall(KSPGetSolution(ksp, &sol));

    510:   /* Get the internal reference to the vector arrays */
    511:   PetscCall(VecGetArrayRead(sol, &x));
    512:   PetscCall(VecGetSize(sol, &N));
    513:   if (err) {
    514:     /* reset the error vector */
    515:     PetscCall(VecSet(err, 0.0));
    516:     /* get array reference */
    517:     PetscCall(VecGetArray(err, &e));
    518:   }

    520:   PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL));

    522:   /* Compute function over the locally owned part of the grid */
    523:   for (moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
    524:     const moab::EntityHandle vhandle = *iter;

    526:     /* get the local DoF numbers to appropriately set the element contribution in the operator */
    527: #ifdef LOCAL_ASSEMBLY
    528:     PetscCall(DMMoabGetFieldDofsLocal(dm, 1, &vhandle, 0, &dof_index));
    529: #else
    530:     PetscCall(DMMoabGetFieldDofs(dm, 1, &vhandle, 0, &dof_index));
    531: #endif

    533:     /* compute the mid-point of the element and use a 1-point lumped quadrature */
    534:     PetscCall(DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos));

    536:     /* compute the discrete L2 error against the exact solution */
    537:     const PetscScalar lerr = (ExactSolution(vpos, user) - x[dof_index]);
    538:     l2err += PetscRealPart(lerr * PetscConj(lerr));
    539:     if (linferr < PetscAbsScalar(lerr)) linferr = PetscAbsScalar(lerr);

    541:     if (err) { /* set the discrete L2 error against the exact solution */
    542:       e[dof_index] = lerr;
    543:     }
    544:   }

    546:   PetscCallMPI(MPIU_Allreduce(&l2err, &global_l2, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD));
    547:   PetscCallMPI(MPIU_Allreduce(&linferr, &global_linf, 1, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD));
    548:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Computed Errors: L_2 = %f, L_inf = %f\n", sqrt(global_l2 / N), global_linf));

    550:   /* Restore vectors */
    551:   PetscCall(VecRestoreArrayRead(sol, &x));
    552:   if (err) PetscCall(VecRestoreArray(err, &e));
    553:   PetscFunctionReturn(PETSC_SUCCESS);
    554: }

    556: PetscErrorCode InitializeOptions(UserContext *user)
    557: {
    558:   const char *bcTypes[2] = {"dirichlet", "neumann"};
    559:   PetscInt    bc;

    561:   PetscFunctionBegin;
    562:   /* set default parameters */
    563:   user->dim       = 2;
    564:   user->problem   = 1;
    565:   user->n         = 2;
    566:   user->nlevels   = 2;
    567:   user->rho       = 0.1;
    568:   user->bounds[0] = user->bounds[2] = user->bounds[4] = 0.0;
    569:   user->bounds[1] = user->bounds[3] = user->bounds[5] = 1.0;
    570:   user->xref                                          = user->bounds[1] / 2;
    571:   user->yref                                          = user->bounds[3] / 2;
    572:   user->nu                                            = 0.05;
    573:   user->usemg                                         = PETSC_FALSE;
    574:   user->io                                            = PETSC_FALSE;
    575:   user->usetri                                        = PETSC_FALSE;
    576:   user->error                                         = PETSC_FALSE;
    577:   bc                                                  = (PetscInt)DIRICHLET;

    579:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "ex35.cxx");
    580:   PetscCall(PetscOptionsInt("-problem", "The type of problem being solved (controls forcing function)", "ex35.cxx", user->problem, &user->problem, NULL));
    581:   PetscCall(PetscOptionsInt("-n", "The elements in each direction", "ex35.cxx", user->n, &user->n, NULL));
    582:   PetscCall(PetscOptionsInt("-levels", "Number of levels in the multigrid hierarchy", "ex35.cxx", user->nlevels, &user->nlevels, NULL));
    583:   PetscCall(PetscOptionsReal("-rho", "The conductivity", "ex35.cxx", user->rho, &user->rho, NULL));
    584:   PetscCall(PetscOptionsReal("-x", "The domain size in x-direction", "ex35.cxx", user->bounds[1], &user->bounds[1], NULL));
    585:   PetscCall(PetscOptionsReal("-y", "The domain size in y-direction", "ex35.cxx", user->bounds[3], &user->bounds[3], NULL));
    586:   PetscCall(PetscOptionsReal("-xref", "The x-coordinate of Gaussian center (for -problem 1)", "ex35.cxx", user->xref, &user->xref, NULL));
    587:   PetscCall(PetscOptionsReal("-yref", "The y-coordinate of Gaussian center (for -problem 1)", "ex35.cxx", user->yref, &user->yref, NULL));
    588:   PetscCall(PetscOptionsReal("-nu", "The width of the Gaussian source (for -problem 1)", "ex35.cxx", user->nu, &user->nu, NULL));
    589:   PetscCall(PetscOptionsBool("-mg", "Use multigrid preconditioner", "ex35.cxx", user->usemg, &user->usemg, NULL));
    590:   PetscCall(PetscOptionsBool("-io", "Write out the solution and mesh data", "ex35.cxx", user->io, &user->io, NULL));
    591:   PetscCall(PetscOptionsBool("-tri", "Use triangles to discretize the domain", "ex35.cxx", user->usetri, &user->usetri, NULL));
    592:   PetscCall(PetscOptionsBool("-error", "Compute the discrete L_2 and L_inf errors of the solution", "ex35.cxx", user->error, &user->error, NULL));
    593:   PetscCall(PetscOptionsEList("-bc", "Type of boundary condition", "ex35.cxx", bcTypes, 2, bcTypes[0], &bc, NULL));
    594:   PetscCall(PetscOptionsString("-file", "The mesh file for the problem", "ex35.cxx", "", user->filename, sizeof(user->filename), &user->use_extfile));
    595:   PetscOptionsEnd();

    597:   if (user->problem < 1 || user->problem > 3) user->problem = 1;
    598:   user->bcType = (BCType)bc;
    599:   user->VPERE  = (user->usetri ? 3 : 4);
    600:   if (user->problem == 3) {
    601:     user->bounds[0] = user->bounds[2] = -0.5;
    602:     user->bounds[1] = user->bounds[3] = 0.5;
    603:     user->bounds[4] = user->bounds[5] = 0.5;
    604:   }
    605:   PetscFunctionReturn(PETSC_SUCCESS);
    606: }

    608: /*TEST

    610:    build:
    611:       requires: moab !complex

    613:    test:
    614:       args: -levels 0 -nu .01 -n 10 -ksp_type cg -pc_type sor -ksp_converged_reason

    616:    test:
    617:       suffix: 2
    618:       nsize: 2
    619:       requires: hdf5
    620:       args: -levels 3 -nu .01 -n 2 -mg -ksp_converged_reason

    622:    test:
    623:       suffix: 3
    624:       nsize: 2
    625:       requires: hdf5
    626:       args: -problem 3 -file data/ex35_mesh.h5m -mg -levels 1 -ksp_converged_reason
    627:       localrunfiles: data

    629: TEST*/