1 static char help[] = "Tests multifield and multicomponent L2 projection.\n"; 2 3 #include <petscdmswarm.h> 4 #include <petscksp.h> 5 #include <petscdmplex.h> 6 #include <petscds.h> 7 8 typedef struct { 9 PetscBool pass; // Don't fail when moments are wrong 10 PetscBool fv; // Use an FV discretization, instead of FE 11 PetscInt Npc; // The number of partices per cell 12 PetscInt field; // The field to project 13 PetscInt Nm; // The number of moments to match 14 PetscReal mtol; // Tolerance for checking moment conservation 15 PetscSimplePointFn *func; // Function used to set particle weights 16 } AppCtx; 17 18 typedef enum { 19 FUNCTION_CONSTANT, 20 FUNCTION_LINEAR, 21 FUNCTION_SIN, 22 FUNCTION_X2_X4, 23 FUNCTION_UNKNOWN, 24 NUM_FUNCTIONS 25 } FunctionType; 26 const char *const FunctionTypes[] = {"constant", "linear", "sin", "x2_x4", "unknown", "FunctionType", "FUNCTION_", NULL}; 27 28 static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 29 { 30 u[0] = 1.0; 31 return PETSC_SUCCESS; 32 } 33 34 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 35 { 36 u[0] = 0.0; 37 for (PetscInt d = 0; d < dim; ++d) u[0] += x[d]; 38 return PETSC_SUCCESS; 39 } 40 41 static PetscErrorCode sinx(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 42 { 43 u[0] = 1; 44 for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSinReal(2. * PETSC_PI * x[d]); 45 return PETSC_SUCCESS; 46 } 47 48 static PetscErrorCode x2_x4(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 49 { 50 u[0] = 1; 51 for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSqr(x[d]) - PetscSqr(PetscSqr(x[d])); 52 return PETSC_SUCCESS; 53 } 54 55 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 56 { 57 FunctionType func = FUNCTION_LINEAR; 58 PetscBool flg; 59 60 PetscFunctionBeginUser; 61 options->pass = PETSC_FALSE; 62 options->fv = PETSC_FALSE; 63 options->Npc = 1; 64 options->field = 0; 65 options->Nm = 1; 66 options->mtol = 100. * PETSC_MACHINE_EPSILON; 67 68 PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX"); 69 PetscCall(PetscOptionsBool("-pass", "Don't fail when moments are wrong", __FILE__, options->pass, &options->pass, NULL)); 70 PetscCall(PetscOptionsBool("-fv", "Use FV instead of FE", __FILE__, options->fv, &options->fv, NULL)); 71 PetscCall(PetscOptionsBoundedInt("-npc", "Number of particles per cell", __FILE__, options->Npc, &options->Npc, NULL, 1)); 72 PetscCall(PetscOptionsBoundedInt("-field", "The field to project", __FILE__, options->field, &options->field, NULL, 0)); 73 PetscCall(PetscOptionsBoundedInt("-moments", "Number of moments to match", __FILE__, options->Nm, &options->Nm, NULL, 0)); 74 PetscCheck(options->Nm < 4, comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot match %" PetscInt_FMT " > 3 moments", options->Nm); 75 PetscCall(PetscOptionsReal("-mtol", "Tolerance for moment checks", "ex2.c", options->mtol, &options->mtol, NULL)); 76 PetscCall(PetscOptionsEnum("-func", "Type of particle weight function", __FILE__, FunctionTypes, (PetscEnum)func, (PetscEnum *)&func, &flg)); 77 switch (func) { 78 case FUNCTION_CONSTANT: 79 options->func = constant; 80 break; 81 case FUNCTION_LINEAR: 82 options->func = linear; 83 break; 84 case FUNCTION_SIN: 85 options->func = sinx; 86 break; 87 case FUNCTION_X2_X4: 88 options->func = x2_x4; 89 break; 90 default: 91 PetscCheck(flg, comm, PETSC_ERR_ARG_WRONG, "Cannnot handle function \"%s\"", FunctionTypes[func]); 92 } 93 PetscOptionsEnd(); 94 95 PetscFunctionReturn(PETSC_SUCCESS); 96 } 97 98 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 99 { 100 PetscFunctionBeginUser; 101 PetscCall(DMCreate(comm, dm)); 102 PetscCall(DMSetType(*dm, DMPLEX)); 103 PetscCall(DMSetFromOptions(*dm)); 104 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107 108 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 109 { 110 PetscFE fe; 111 PetscFV fv; 112 DMPolytopeType ct; 113 PetscInt dim, cStart; 114 115 PetscFunctionBeginUser; 116 PetscCall(DMGetDimension(dm, &dim)); 117 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 118 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 119 if (user->fv) { 120 PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv)); 121 PetscCall(PetscObjectSetName((PetscObject)fv, "fv")); 122 PetscCall(PetscFVSetNumComponents(fv, 1)); 123 PetscCall(PetscFVSetSpatialDimension(fv, dim)); 124 PetscCall(PetscFVCreateDualSpace(fv, ct)); 125 PetscCall(PetscFVSetFromOptions(fv)); 126 PetscCall(DMAddField(dm, NULL, (PetscObject)fv)); 127 PetscCall(PetscFVDestroy(&fv)); 128 PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv)); 129 PetscCall(PetscObjectSetName((PetscObject)fv, "fv2")); 130 PetscCall(PetscFVSetNumComponents(fv, 2)); 131 PetscCall(PetscFVSetSpatialDimension(fv, dim)); 132 PetscCall(PetscFVCreateDualSpace(fv, ct)); 133 PetscCall(PetscFVSetFromOptions(fv)); 134 PetscCall(DMAddField(dm, NULL, (PetscObject)fv)); 135 PetscCall(PetscFVDestroy(&fv)); 136 } else { 137 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, -1, &fe)); 138 PetscCall(PetscObjectSetName((PetscObject)fe, "fe")); 139 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 140 PetscCall(PetscFEDestroy(&fe)); 141 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 2, ct, NULL, -1, &fe)); 142 PetscCall(PetscObjectSetName((PetscObject)fe, "fe2")); 143 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 144 PetscCall(PetscFEDestroy(&fe)); 145 } 146 PetscCall(DMCreateDS(dm)); 147 if (user->fv) { 148 DMLabel label; 149 PetscInt values[1] = {1}; 150 151 PetscCall(DMCreateLabel(dm, "ghost")); 152 PetscCall(DMGetLabel(dm, "marker", &label)); 153 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "dummy", label, 1, values, 0, 0, NULL, NULL, NULL, NULL, NULL)); 154 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "dummy", label, 1, values, 1, 0, NULL, NULL, NULL, NULL, NULL)); 155 } 156 PetscFunctionReturn(PETSC_SUCCESS); 157 } 158 159 static PetscErrorCode CreateSwarm(DM dm, DM *sw, AppCtx *user) 160 { 161 PetscScalar *coords, *wvals, *xvals; 162 PetscInt Npc = user->Npc, dim, Np; 163 164 PetscFunctionBeginUser; 165 PetscCall(DMGetDimension(dm, &dim)); 166 167 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw)); 168 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles")); 169 PetscCall(DMSetType(*sw, DMSWARM)); 170 PetscCall(DMSetDimension(*sw, dim)); 171 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 172 PetscCall(DMSwarmSetCellDM(*sw, dm)); 173 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR)); 174 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "x_q", 2, PETSC_SCALAR)); 175 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 176 PetscCall(DMSwarmInsertPointsUsingCellDM(*sw, DMSWARMPIC_LAYOUT_GAUSS, Npc)); 177 PetscCall(DMSetFromOptions(*sw)); 178 179 PetscCall(DMSwarmGetLocalSize(*sw, &Np)); 180 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 181 PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&wvals)); 182 PetscCall(DMSwarmGetField(*sw, "x_q", NULL, NULL, (void **)&xvals)); 183 for (PetscInt p = 0; p < Np; ++p) { 184 PetscCall(user->func(dim, 0., &coords[p * dim], 1, &wvals[p], user)); 185 for (PetscInt c = 0; c < 2; ++c) PetscCall(user->func(dim, 0., &coords[p * dim], 1, &xvals[p * 2 + c], user)); 186 } 187 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 188 PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&wvals)); 189 PetscCall(DMSwarmRestoreField(*sw, "x_q", NULL, NULL, (void **)&xvals)); 190 191 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view")); 192 PetscFunctionReturn(PETSC_SUCCESS); 193 } 194 195 static PetscErrorCode computeParticleMoments(DM sw, Vec u, PetscReal moments[3], AppCtx *user) 196 { 197 DM dm; 198 const PetscReal *coords; 199 const PetscScalar *w; 200 PetscReal mom[3] = {0.0, 0.0, 0.0}; 201 PetscInt dim, cStart, cEnd, Nc; 202 203 PetscFunctionBeginUser; 204 PetscCall(DMGetDimension(sw, &dim)); 205 PetscCall(DMSwarmGetCellDM(sw, &dm)); 206 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 207 PetscCall(DMSwarmSortGetAccess(sw)); 208 PetscCall(DMSwarmGetFieldInfo(sw, user->field ? "x_q" : "w_q", &Nc, NULL)); 209 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 210 PetscCall(VecGetArrayRead(u, &w)); 211 for (PetscInt cell = cStart; cell < cEnd; ++cell) { 212 PetscInt *pidx; 213 PetscInt Np; 214 215 PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx)); 216 for (PetscInt p = 0; p < Np; ++p) { 217 const PetscInt idx = pidx[p]; 218 const PetscReal *x = &coords[idx * dim]; 219 220 for (PetscInt c = 0; c < Nc; ++c) { 221 mom[0] += PetscRealPart(w[idx * Nc + c]); 222 mom[1] += PetscRealPart(w[idx * Nc + c]) * x[0]; 223 for (PetscInt d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx * Nc + c]) * PetscSqr(x[d]); 224 } 225 } 226 PetscCall(PetscFree(pidx)); 227 } 228 PetscCall(VecRestoreArrayRead(u, &w)); 229 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 230 PetscCall(DMSwarmSortRestoreAccess(sw)); 231 PetscCall(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw))); 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 static void f0_1(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 236 { 237 const PetscInt Nc = uOff[1] - uOff[0]; 238 239 for (PetscInt c = 0; c < Nc; ++c) f0[0] += u[c]; 240 } 241 242 static void f0_x(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 243 { 244 const PetscInt Nc = uOff[1] - uOff[0]; 245 246 for (PetscInt c = 0; c < Nc; ++c) f0[0] += x[0] * u[c]; 247 } 248 249 static void f0_r2(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 250 { 251 const PetscInt Nc = uOff[1] - uOff[0]; 252 253 for (PetscInt c = 0; c < Nc; ++c) 254 for (PetscInt d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[c]; 255 } 256 257 static PetscErrorCode computeFieldMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user) 258 { 259 PetscDS ds; 260 PetscScalar mom; 261 262 PetscFunctionBeginUser; 263 PetscCall(DMGetDS(dm, &ds)); 264 PetscCall(PetscDSSetObjective(ds, 0, &f0_1)); 265 mom = 0.; 266 PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user)); 267 moments[0] = PetscRealPart(mom); 268 PetscCall(PetscDSSetObjective(ds, 0, &f0_x)); 269 mom = 0.; 270 PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user)); 271 moments[1] = PetscRealPart(mom); 272 PetscCall(PetscDSSetObjective(ds, 0, &f0_r2)); 273 mom = 0.; 274 PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user)); 275 moments[2] = PetscRealPart(mom); 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 static PetscErrorCode TestParticlesToField(DM sw, DM dm, Vec fhat, AppCtx *user) 280 { 281 const char *fieldnames[1] = {user->field ? "x_q" : "w_q"}; 282 Vec fields[1] = {fhat}, f; 283 PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f 284 PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f 285 286 PetscFunctionBeginUser; 287 PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_FORWARD)); 288 289 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fieldnames[0], &f)); 290 PetscCall(computeParticleMoments(sw, f, pmoments, user)); 291 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fieldnames[0], &f)); 292 PetscCall(computeFieldMoments(dm, fhat, fmoments, user)); 293 PetscCall(VecViewFromOptions(fhat, NULL, "-fhat_view")); 294 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2])); 295 for (PetscInt m = 0; m < user->Nm; ++m) { 296 if (user->pass) { 297 if (PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) > user->mtol) { 298 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "p projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", pmoments[0], pmoments[1], pmoments[2])); 299 } 300 } else { 301 PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->mtol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]), 302 user->mtol); 303 } 304 } 305 PetscFunctionReturn(PETSC_SUCCESS); 306 } 307 308 static PetscErrorCode TestFieldToParticles(DM sw, DM dm, Vec fhat, AppCtx *user) 309 { 310 const char *fieldnames[1] = {user->field ? "x_q" : "w_q"}; 311 Vec fields[1] = {fhat}, f; 312 PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f 313 PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f 314 315 PetscFunctionBeginUser; 316 PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_REVERSE)); 317 318 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fieldnames[0], &f)); 319 PetscCall(computeParticleMoments(sw, f, pmoments, user)); 320 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fieldnames[0], &f)); 321 PetscCall(computeFieldMoments(dm, fhat, fmoments, user)); 322 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2])); 323 for (PetscInt m = 0; m < user->Nm; ++m) { 324 if (user->pass) { 325 if (PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) > user->mtol) { 326 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "p projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", pmoments[0], pmoments[1], pmoments[2])); 327 } 328 } else { 329 PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->mtol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]), 330 user->mtol); 331 } 332 } 333 PetscFunctionReturn(PETSC_SUCCESS); 334 } 335 336 int main(int argc, char *argv[]) 337 { 338 DM dm, subdm, sw; 339 Vec fhat; 340 IS subis; 341 AppCtx user; 342 343 PetscFunctionBeginUser; 344 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 345 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 346 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &user)); 347 PetscCall(CreateDiscretization(dm, &user)); 348 PetscCall(CreateSwarm(dm, &sw, &user)); 349 350 PetscCall(DMCreateSubDM(dm, 1, &user.field, &subis, &subdm)); 351 352 PetscCall(DMGetGlobalVector(subdm, &fhat)); 353 PetscCall(PetscObjectSetName((PetscObject)fhat, "FEM f")); 354 PetscCall(TestParticlesToField(sw, subdm, fhat, &user)); 355 PetscCall(TestFieldToParticles(sw, subdm, fhat, &user)); 356 PetscCall(DMRestoreGlobalVector(subdm, &fhat)); 357 358 PetscCall(ISDestroy(&subis)); 359 PetscCall(DMDestroy(&subdm)); 360 PetscCall(DMDestroy(&dm)); 361 PetscCall(DMDestroy(&sw)); 362 PetscCall(PetscFinalize()); 363 return 0; 364 } 365 366 /*TEST 367 368 # Swarm does not handle complex or quad 369 build: 370 requires: !complex double 371 372 testset: 373 requires: triangle 374 args: -dm_refine 1 -petscspace_degree 2 -moments 3 \ 375 -ptof_pc_type lu \ 376 -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none 377 378 test: 379 suffix: tri_fe_f0 380 args: -field 0 381 382 test: 383 suffix: tri_fe_f1 384 args: -field 1 385 386 test: 387 suffix: quad_fe_f0 388 args: -dm_plex_simplex 0 -field 0 389 390 test: 391 suffix: quad_fe_f1 392 args: -dm_plex_simplex 0 -field 1 393 394 testset: 395 requires: triangle 396 args: -dm_refine 1 -moments 1 -fv \ 397 -ptof_pc_type lu \ 398 -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none 399 400 test: 401 suffix: tri_fv_f0 402 args: -field 0 403 404 test: 405 suffix: tri_fv_f1 406 args: -field 1 407 408 test: 409 suffix: quad_fv_f0 410 args: -dm_plex_simplex 0 -field 0 411 412 test: 413 suffix: quad_fv_f1 414 args: -dm_plex_simplex 0 -field 1 415 416 TEST*/ 417