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