1 static char help[] = "Tests projection with DMSwarm using general particle shapes.\n"; 2 3 #include <petsc/private/dmswarmimpl.h> 4 #include <petsc/private/petscfeimpl.h> 5 6 #include <petscdmplex.h> 7 #include <petscds.h> 8 #include <petscksp.h> 9 10 typedef struct { 11 PetscInt dummy; 12 } AppCtx; 13 14 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) { 15 PetscFunctionBeginUser; 16 PetscCall(DMCreate(comm, dm)); 17 PetscCall(DMSetType(*dm, DMPLEX)); 18 PetscCall(DMSetFromOptions(*dm)); 19 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 20 PetscFunctionReturn(0); 21 } 22 23 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 24 PetscInt d; 25 u[0] = 0.0; 26 for (d = 0; d < dim; ++d) u[0] += x[d]; 27 return 0; 28 } 29 30 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[]) { 31 g0[0] = 1.0; 32 } 33 34 static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user) { 35 PetscDS prob; 36 PetscFE fe; 37 PetscQuadrature quad; 38 PetscScalar *vals; 39 PetscReal *v0, *J, *invJ, detJ, *coords, *xi0; 40 PetscInt *cellid; 41 const PetscReal *qpoints; 42 PetscInt Ncell, c, Nq, q, dim; 43 PetscBool simplex; 44 45 PetscFunctionBeginUser; 46 PetscCall(DMGetDimension(dm, &dim)); 47 PetscCall(DMPlexIsSimplex(dm, &simplex)); 48 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe)); 49 PetscCall(DMGetDS(dm, &prob)); 50 PetscCall(PetscDSSetDiscretization(prob, 0, (PetscObject)fe)); 51 PetscCall(PetscFEDestroy(&fe)); 52 PetscCall(PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL)); 53 PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Ncell)); 54 PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe)); 55 PetscCall(PetscFEGetQuadrature(fe, &quad)); 56 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL)); 57 58 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw)); 59 PetscCall(DMSetType(*sw, DMSWARM)); 60 PetscCall(DMSetDimension(*sw, dim)); 61 62 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 63 PetscCall(DMSwarmSetCellDM(*sw, dm)); 64 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_SCALAR)); 65 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 66 PetscCall(DMSwarmSetLocalSizes(*sw, Ncell * Nq, 0)); 67 PetscCall(DMSetFromOptions(*sw)); 68 69 PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ)); 70 for (c = 0; c < dim; c++) xi0[c] = -1.; 71 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 72 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid)); 73 PetscCall(DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **)&vals)); 74 for (c = 0; c < Ncell; ++c) { 75 for (q = 0; q < Nq; ++q) { 76 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 77 cellid[c * Nq + q] = c; 78 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]); 79 linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], NULL); 80 } 81 } 82 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 83 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid)); 84 PetscCall(DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **)&vals)); 85 PetscCall(PetscFree4(xi0, v0, J, invJ)); 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user) { 90 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 91 KSP ksp; 92 Mat mass; 93 Vec u, rhs, uproj; 94 PetscReal error; 95 96 PetscFunctionBeginUser; 97 funcs[0] = linear; 98 99 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u)); 100 PetscCall(VecViewFromOptions(u, NULL, "-f_view")); 101 PetscCall(DMGetGlobalVector(dm, &rhs)); 102 PetscCall(DMCreateMassMatrix(sw, dm, &mass)); 103 PetscCall(MatMult(mass, u, rhs)); 104 PetscCall(MatDestroy(&mass)); 105 PetscCall(VecDestroy(&u)); 106 107 PetscCall(DMGetGlobalVector(dm, &uproj)); 108 PetscCall(DMCreateMatrix(dm, &mass)); 109 PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user)); 110 PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view")); 111 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 112 PetscCall(KSPSetOperators(ksp, mass, mass)); 113 PetscCall(KSPSetFromOptions(ksp)); 114 PetscCall(KSPSolve(ksp, rhs, uproj)); 115 PetscCall(KSPDestroy(&ksp)); 116 PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection")); 117 PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view")); 118 PetscCall(DMComputeL2Diff(dm, 0.0, funcs, NULL, uproj, &error)); 119 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error)); 120 PetscCall(MatDestroy(&mass)); 121 PetscCall(DMRestoreGlobalVector(dm, &rhs)); 122 PetscCall(DMRestoreGlobalVector(dm, &uproj)); 123 PetscFunctionReturn(0); 124 } 125 126 int main(int argc, char *argv[]) { 127 MPI_Comm comm; 128 DM dm, sw; 129 AppCtx user; 130 131 PetscFunctionBeginUser; 132 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 133 comm = PETSC_COMM_WORLD; 134 PetscCall(CreateMesh(comm, &dm, &user)); 135 PetscCall(CreateParticles(dm, &sw, &user)); 136 PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh")); 137 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 138 PetscCall(DMViewFromOptions(sw, NULL, "-sw_view")); 139 PetscCall(TestL2Projection(dm, sw, &user)); 140 PetscCall(DMDestroy(&dm)); 141 PetscCall(DMDestroy(&sw)); 142 PetscFinalize(); 143 return 0; 144 } 145 146 /*TEST 147 148 test: 149 suffix: proj_0 150 requires: pragmatic 151 TODO: broken 152 args: -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 153 test: 154 suffix: proj_1 155 requires: pragmatic 156 TODO: broken 157 args: -dm_plex_simplex 0 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 158 test: 159 suffix: proj_2 160 requires: pragmatic 161 TODO: broken 162 args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 163 test: 164 suffix: proj_3 165 requires: pragmatic 166 TODO: broken 167 args: -dm_plex_simplex 0 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 168 169 TEST*/ 170