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