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, "Cannot handle function \"%s\"", FunctionTypes[func]); 92 } 93 PetscOptionsEnd(); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 98 { 99 PetscFunctionBeginUser; 100 PetscCall(DMCreate(comm, dm)); 101 PetscCall(DMSetType(*dm, DMPLEX)); 102 PetscCall(DMSetFromOptions(*dm)); 103 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 104 PetscFunctionReturn(PETSC_SUCCESS); 105 } 106 107 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 108 { 109 PetscFE fe; 110 PetscFV fv; 111 DMPolytopeType ct; 112 PetscInt dim, cStart; 113 114 PetscFunctionBeginUser; 115 PetscCall(DMGetDimension(dm, &dim)); 116 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 117 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 118 if (user->fv) { 119 PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv)); 120 PetscCall(PetscObjectSetName((PetscObject)fv, "fv")); 121 PetscCall(PetscFVSetNumComponents(fv, 1)); 122 PetscCall(PetscFVSetSpatialDimension(fv, dim)); 123 PetscCall(PetscFVCreateDualSpace(fv, ct)); 124 PetscCall(PetscFVSetFromOptions(fv)); 125 PetscCall(DMAddField(dm, NULL, (PetscObject)fv)); 126 PetscCall(PetscFVDestroy(&fv)); 127 PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv)); 128 PetscCall(PetscObjectSetName((PetscObject)fv, "fv2")); 129 PetscCall(PetscFVSetNumComponents(fv, 2)); 130 PetscCall(PetscFVSetSpatialDimension(fv, dim)); 131 PetscCall(PetscFVCreateDualSpace(fv, ct)); 132 PetscCall(PetscFVSetFromOptions(fv)); 133 PetscCall(DMAddField(dm, NULL, (PetscObject)fv)); 134 PetscCall(PetscFVDestroy(&fv)); 135 } else { 136 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, -1, &fe)); 137 PetscCall(PetscObjectSetName((PetscObject)fe, "fe")); 138 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 139 PetscCall(PetscFEDestroy(&fe)); 140 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 2, ct, NULL, -1, &fe)); 141 PetscCall(PetscObjectSetName((PetscObject)fe, "fe2")); 142 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 143 PetscCall(PetscFEDestroy(&fe)); 144 } 145 PetscCall(DMCreateDS(dm)); 146 if (user->fv) { 147 DMLabel label; 148 PetscInt values[1] = {1}; 149 150 PetscCall(DMCreateLabel(dm, "ghost")); 151 PetscCall(DMGetLabel(dm, "marker", &label)); 152 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "dummy", label, 1, values, 0, 0, NULL, NULL, NULL, NULL, NULL)); 153 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "dummy", label, 1, values, 1, 0, NULL, NULL, NULL, NULL, NULL)); 154 } 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 static PetscErrorCode CreateSwarm(DM dm, DM *sw, AppCtx *user) 159 { 160 PetscScalar *coords, *wvals, *xvals; 161 PetscInt Npc = user->Npc, dim, Np; 162 163 PetscFunctionBeginUser; 164 PetscCall(DMGetDimension(dm, &dim)); 165 166 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw)); 167 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles")); 168 PetscCall(DMSetType(*sw, DMSWARM)); 169 PetscCall(DMSetDimension(*sw, dim)); 170 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 171 PetscCall(DMSwarmSetCellDM(*sw, dm)); 172 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR)); 173 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "x_q", 2, PETSC_SCALAR)); 174 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 175 PetscCall(DMSwarmInsertPointsUsingCellDM(*sw, DMSWARMPIC_LAYOUT_GAUSS, Npc)); 176 PetscCall(DMSetFromOptions(*sw)); 177 178 PetscCall(DMSwarmGetLocalSize(*sw, &Np)); 179 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 180 PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&wvals)); 181 PetscCall(DMSwarmGetField(*sw, "x_q", NULL, NULL, (void **)&xvals)); 182 for (PetscInt p = 0; p < Np; ++p) { 183 PetscCall(user->func(dim, 0., &coords[p * dim], 1, &wvals[p], user)); 184 for (PetscInt c = 0; c < 2; ++c) PetscCall(user->func(dim, 0., &coords[p * dim], 1, &xvals[p * 2 + c], user)); 185 } 186 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 187 PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&wvals)); 188 PetscCall(DMSwarmRestoreField(*sw, "x_q", NULL, NULL, (void **)&xvals)); 189 190 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view")); 191 PetscFunctionReturn(PETSC_SUCCESS); 192 } 193 194 static PetscErrorCode computeParticleMoments(DM sw, Vec u, PetscReal moments[3], AppCtx *user) 195 { 196 DM dm; 197 const PetscReal *coords; 198 const PetscScalar *w; 199 PetscReal mom[3] = {0.0, 0.0, 0.0}; 200 PetscInt dim, cStart, cEnd, Nc; 201 202 PetscFunctionBeginUser; 203 PetscCall(DMGetDimension(sw, &dim)); 204 PetscCall(DMSwarmGetCellDM(sw, &dm)); 205 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 206 PetscCall(DMSwarmSortGetAccess(sw)); 207 PetscCall(DMSwarmGetFieldInfo(sw, user->field ? "x_q" : "w_q", &Nc, NULL)); 208 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 209 PetscCall(VecGetArrayRead(u, &w)); 210 for (PetscInt cell = cStart; cell < cEnd; ++cell) { 211 PetscInt *pidx; 212 PetscInt Np; 213 214 PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx)); 215 for (PetscInt p = 0; p < Np; ++p) { 216 const PetscInt idx = pidx[p]; 217 const PetscReal *x = &coords[idx * dim]; 218 219 for (PetscInt c = 0; c < Nc; ++c) { 220 mom[0] += PetscRealPart(w[idx * Nc + c]); 221 mom[1] += PetscRealPart(w[idx * Nc + c]) * x[0]; 222 for (PetscInt d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx * Nc + c]) * PetscSqr(x[d]); 223 } 224 } 225 PetscCall(PetscFree(pidx)); 226 } 227 PetscCall(VecRestoreArrayRead(u, &w)); 228 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 229 PetscCall(DMSwarmSortRestoreAccess(sw)); 230 PetscCall(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw))); 231 PetscFunctionReturn(PETSC_SUCCESS); 232 } 233 234 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[]) 235 { 236 const PetscInt Nc = uOff[1] - uOff[0]; 237 238 for (PetscInt c = 0; c < Nc; ++c) f0[0] += u[c]; 239 } 240 241 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[]) 242 { 243 const PetscInt Nc = uOff[1] - uOff[0]; 244 245 for (PetscInt c = 0; c < Nc; ++c) f0[0] += x[0] * u[c]; 246 } 247 248 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[]) 249 { 250 const PetscInt Nc = uOff[1] - uOff[0]; 251 252 for (PetscInt c = 0; c < Nc; ++c) 253 for (PetscInt d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[c]; 254 } 255 256 static PetscErrorCode computeFieldMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user) 257 { 258 PetscDS ds; 259 PetscScalar mom; 260 261 PetscFunctionBeginUser; 262 PetscCall(DMGetDS(dm, &ds)); 263 PetscCall(PetscDSSetObjective(ds, 0, &f0_1)); 264 mom = 0.; 265 PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user)); 266 moments[0] = PetscRealPart(mom); 267 PetscCall(PetscDSSetObjective(ds, 0, &f0_x)); 268 mom = 0.; 269 PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user)); 270 moments[1] = PetscRealPart(mom); 271 PetscCall(PetscDSSetObjective(ds, 0, &f0_r2)); 272 mom = 0.; 273 PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user)); 274 moments[2] = PetscRealPart(mom); 275 PetscFunctionReturn(PETSC_SUCCESS); 276 } 277 278 static PetscErrorCode TestParticlesToField(DM sw, DM dm, Vec fhat, AppCtx *user) 279 { 280 const char *fieldnames[1] = {user->field ? "x_q" : "w_q"}; 281 Vec fields[1] = {fhat}, f; 282 PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f 283 PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f 284 285 PetscFunctionBeginUser; 286 PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_FORWARD)); 287 288 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fieldnames[0], &f)); 289 PetscCall(computeParticleMoments(sw, f, pmoments, user)); 290 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fieldnames[0], &f)); 291 PetscCall(computeFieldMoments(dm, fhat, fmoments, user)); 292 PetscCall(VecViewFromOptions(fhat, NULL, "-fhat_view")); 293 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2])); 294 for (PetscInt m = 0; m < user->Nm; ++m) { 295 if (user->pass) { 296 if (PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) > user->mtol) { 297 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "p projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", pmoments[0], pmoments[1], pmoments[2])); 298 } 299 } else { 300 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]), 301 user->mtol); 302 } 303 } 304 PetscFunctionReturn(PETSC_SUCCESS); 305 } 306 307 static PetscErrorCode TestFieldToParticles(DM sw, DM dm, Vec fhat, AppCtx *user) 308 { 309 const char *fieldnames[1] = {user->field ? "x_q" : "w_q"}; 310 Vec fields[1] = {fhat}, f; 311 PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f 312 PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f 313 314 PetscFunctionBeginUser; 315 PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_REVERSE)); 316 317 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fieldnames[0], &f)); 318 PetscCall(computeParticleMoments(sw, f, pmoments, user)); 319 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fieldnames[0], &f)); 320 PetscCall(computeFieldMoments(dm, fhat, fmoments, user)); 321 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2])); 322 for (PetscInt m = 0; m < user->Nm; ++m) { 323 if (user->pass) { 324 if (PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) > user->mtol) { 325 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "p projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", pmoments[0], pmoments[1], pmoments[2])); 326 } 327 } else { 328 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]), 329 user->mtol); 330 } 331 } 332 PetscFunctionReturn(PETSC_SUCCESS); 333 } 334 335 int main(int argc, char *argv[]) 336 { 337 DM dm, subdm, sw; 338 Vec fhat; 339 IS subis; 340 AppCtx user; 341 342 PetscFunctionBeginUser; 343 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 344 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 345 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &user)); 346 PetscCall(CreateDiscretization(dm, &user)); 347 PetscCall(CreateSwarm(dm, &sw, &user)); 348 349 PetscCall(DMCreateSubDM(dm, 1, &user.field, &subis, &subdm)); 350 351 PetscCall(DMGetGlobalVector(subdm, &fhat)); 352 PetscCall(PetscObjectSetName((PetscObject)fhat, "FEM f")); 353 PetscCall(TestParticlesToField(sw, subdm, fhat, &user)); 354 PetscCall(TestFieldToParticles(sw, subdm, fhat, &user)); 355 PetscCall(DMRestoreGlobalVector(subdm, &fhat)); 356 357 PetscCall(ISDestroy(&subis)); 358 PetscCall(DMDestroy(&subdm)); 359 PetscCall(DMDestroy(&dm)); 360 PetscCall(DMDestroy(&sw)); 361 PetscCall(PetscFinalize()); 362 return 0; 363 } 364 365 /*TEST 366 367 # Swarm does not handle complex or quad 368 build: 369 requires: !complex double 370 371 testset: 372 requires: triangle 373 args: -dm_refine 1 -petscspace_degree 2 -moments 3 \ 374 -ptof_pc_type lu \ 375 -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none 376 377 test: 378 suffix: tri_fe_f0 379 args: -field 0 380 381 test: 382 suffix: tri_fe_f1 383 args: -field 1 384 385 test: 386 suffix: quad_fe_f0 387 args: -dm_plex_simplex 0 -field 0 388 389 test: 390 suffix: quad_fe_f1 391 args: -dm_plex_simplex 0 -field 1 392 393 testset: 394 requires: triangle 395 args: -dm_refine 1 -moments 1 -fv \ 396 -ptof_pc_type lu \ 397 -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none 398 399 test: 400 suffix: tri_fv_f0 401 args: -field 0 402 403 test: 404 suffix: tri_fv_f1 405 args: -field 1 406 407 test: 408 suffix: quad_fv_f0 409 args: -dm_plex_simplex 0 -field 0 410 411 test: 412 suffix: quad_fv_f1 413 args: -dm_plex_simplex 0 -field 1 414 415 TEST*/ 416