1 static char help[] = "Tests projection with DMSwarm using general particle shapes.\n"; 2 3 #include <petsc/private/petscfeimpl.h> 4 #include <petscdmswarm.h> 5 #include <petscdmplex.h> 6 #include <petscds.h> 7 #include <petscksp.h> 8 9 typedef struct { 10 PetscReal L[3]; /* Dimensions of the mesh bounding box */ 11 } AppCtx; 12 13 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 14 { 15 PetscReal low[3], high[3]; 16 PetscInt dim; 17 18 PetscFunctionBeginUser; 19 PetscCall(DMCreate(comm, dm)); 20 PetscCall(DMSetType(*dm, DMPLEX)); 21 PetscCall(DMSetFromOptions(*dm)); 22 PetscCall(DMGetBoundingBox(*dm, low, high)); 23 PetscCall(DMGetDimension(*dm, &dim)); 24 for (PetscInt d = 0; d < dim; ++d) user->L[d] = high[d] - low[d]; 25 PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh")); 26 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 27 PetscFunctionReturn(PETSC_SUCCESS); 28 } 29 30 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx) 31 { 32 AppCtx *ctx = (AppCtx *)a_ctx; 33 PetscInt d; 34 u[0] = 0.0; 35 for (d = 0; d < dim; ++d) u[0] += x[d] / ctx->L[d]; 36 return PETSC_SUCCESS; 37 } 38 39 static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 40 { 41 g0[0] = 1.0; 42 } 43 44 static PetscErrorCode CreateFEM(DM dm, AppCtx *user) 45 { 46 PetscFE fe; 47 PetscDS ds; 48 PetscBool simplex; 49 PetscInt dim; 50 51 PetscFunctionBeginUser; 52 PetscCall(DMGetDimension(dm, &dim)); 53 PetscCall(DMPlexIsSimplex(dm, &simplex)); 54 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe)); 55 PetscCall(PetscObjectSetName((PetscObject)fe, "fe")); 56 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 57 PetscCall(DMCreateDS(dm)); 58 PetscCall(PetscFEDestroy(&fe)); 59 /* Setup to form mass matrix */ 60 PetscCall(DMGetDS(dm, &ds)); 61 PetscCall(PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL)); 62 PetscFunctionReturn(PETSC_SUCCESS); 63 } 64 65 static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user) 66 { 67 PetscDS prob; 68 PetscFE fe; 69 PetscQuadrature quad; 70 PetscScalar *vals; 71 PetscReal *v0, *J, *invJ, detJ, *coords, *xi0; 72 PetscInt *cellid; 73 const PetscReal *qpoints, *qweights; 74 PetscInt cStart, cEnd, c, Nq, q, dim; 75 76 PetscFunctionBeginUser; 77 PetscCall(DMGetDimension(dm, &dim)); 78 PetscCall(DMGetDS(dm, &prob)); 79 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 80 PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe)); 81 PetscCall(PetscFEGetQuadrature(fe, &quad)); 82 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, &qweights)); 83 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw)); 84 PetscCall(DMSetType(*sw, DMSWARM)); 85 PetscCall(DMSetDimension(*sw, dim)); 86 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 87 PetscCall(DMSwarmSetCellDM(*sw, dm)); 88 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_REAL)); 89 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 90 PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Nq, 0)); 91 PetscCall(DMSetFromOptions(*sw)); 92 PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ)); 93 for (c = 0; c < dim; c++) xi0[c] = -1.; 94 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 95 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid)); 96 PetscCall(DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **)&vals)); 97 for (c = cStart; c < cEnd; ++c) { 98 for (q = 0; q < Nq; ++q) { 99 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 100 cellid[c * Nq + q] = c; 101 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]); 102 PetscCall(linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], (void *)user)); 103 vals[c * Nq + q] *= qweights[q] * detJ; 104 } 105 } 106 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 107 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid)); 108 PetscCall(DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **)&vals)); 109 PetscCall(PetscFree4(xi0, v0, J, invJ)); 110 PetscCall(DMSwarmMigrate(*sw, PETSC_FALSE)); 111 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles")); 112 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view")); 113 PetscFunctionReturn(PETSC_SUCCESS); 114 } 115 116 static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user) 117 { 118 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 119 KSP ksp; 120 Mat mass; 121 Vec u, rhs, uproj; 122 void *ctxs[1]; 123 PetscReal error; 124 125 PetscFunctionBeginUser; 126 ctxs[0] = (void *)user; 127 funcs[0] = linear; 128 129 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u)); 130 PetscCall(VecViewFromOptions(u, NULL, "-f_view")); 131 PetscCall(DMGetGlobalVector(dm, &rhs)); 132 PetscCall(DMCreateMassMatrix(sw, dm, &mass)); 133 PetscCall(MatMultTranspose(mass, u, rhs)); 134 PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view")); 135 PetscCall(MatDestroy(&mass)); 136 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "f_q", &u)); 137 PetscCall(DMGetGlobalVector(dm, &uproj)); 138 PetscCall(DMCreateMatrix(dm, &mass)); 139 PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user)); 140 PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view")); 141 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 142 PetscCall(KSPSetOperators(ksp, mass, mass)); 143 PetscCall(KSPSetFromOptions(ksp)); 144 PetscCall(KSPSolve(ksp, rhs, uproj)); 145 PetscCall(KSPDestroy(&ksp)); 146 PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection")); 147 PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view")); 148 PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, uproj, &error)); 149 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error)); 150 PetscCall(MatDestroy(&mass)); 151 PetscCall(DMRestoreGlobalVector(dm, &rhs)); 152 PetscCall(DMRestoreGlobalVector(dm, &uproj)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 int main(int argc, char *argv[]) 157 { 158 MPI_Comm comm; 159 DM dm, sw; 160 AppCtx user; 161 162 PetscFunctionBeginUser; 163 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 164 comm = PETSC_COMM_WORLD; 165 PetscCall(CreateMesh(comm, &dm, &user)); 166 PetscCall(CreateFEM(dm, &user)); 167 PetscCall(CreateParticles(dm, &sw, &user)); 168 PetscCall(TestL2Projection(dm, sw, &user)); 169 PetscCall(DMDestroy(&dm)); 170 PetscCall(DMDestroy(&sw)); 171 PetscCall(PetscFinalize()); 172 return 0; 173 } 174 175 /*TEST 176 177 build: 178 requires: double !complex 179 test: 180 suffix: proj_0 181 requires: triangle 182 args: -dm_plex_box_faces 2,2\ 183 -dm_view -sw_view\ 184 -petscspace_degree 1 -petscfe_default_quadrature_order 1\ 185 -pc_type lu 186 test: 187 suffix: proj_1 188 requires: triangle 189 args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2\ 190 -dm_view -sw_view\ 191 -petscspace_degree 1 -petscfe_default_quadrature_order 1\ 192 -pc_type lu 193 test: 194 suffix: proj_2 195 requires: ctetgen 196 args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\ 197 -dm_view -sw_view\ 198 -petscspace_degree 1 -petscfe_default_quadrature_order 1\ 199 -pc_type lu 200 test: 201 suffix: proj_3 202 requires: ctetgen 203 args: -dm_plex_dim 3 -dm_plex_simplex 0\ 204 -dm_plex_box_faces 2,2,2\ 205 -dm_view -sw_view\ 206 -petscspace_degree 1 -petscfe_default_quadrature_order 1\ 207 -pc_type lu 208 209 TEST*/ 210