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 CHKERRQ(DMCreate(comm, dm)); 18 CHKERRQ(DMSetType(*dm, DMPLEX)); 19 CHKERRQ(DMSetFromOptions(*dm)); 20 CHKERRQ(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 CHKERRQ(DMGetDimension(dm, &dim)); 54 CHKERRQ(DMPlexIsSimplex(dm, &simplex)); 55 CHKERRQ(PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, NULL, -1, &fe)); 56 CHKERRQ(DMGetDS(dm, &prob)); 57 CHKERRQ(PetscDSSetDiscretization(prob, 0, (PetscObject) fe)); 58 CHKERRQ(PetscFEDestroy(&fe)); 59 CHKERRQ(PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL)); 60 CHKERRQ(DMPlexGetHeightStratum(dm, 0, NULL, &Ncell)); 61 CHKERRQ(PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe)); 62 CHKERRQ(PetscFEGetQuadrature(fe, &quad)); 63 CHKERRQ(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL)); 64 65 CHKERRQ(DMCreate(PetscObjectComm((PetscObject) dm), sw)); 66 CHKERRQ(DMSetType(*sw, DMSWARM)); 67 CHKERRQ(DMSetDimension(*sw, dim)); 68 69 CHKERRQ(DMSwarmSetType(*sw, DMSWARM_PIC)); 70 CHKERRQ(DMSwarmSetCellDM(*sw, dm)); 71 CHKERRQ(DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_SCALAR)); 72 CHKERRQ(DMSwarmFinalizeFieldRegister(*sw)); 73 CHKERRQ(DMSwarmSetLocalSizes(*sw, Ncell * Nq, 0)); 74 CHKERRQ(DMSetFromOptions(*sw)); 75 76 CHKERRQ(PetscMalloc4(dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ)); 77 for (c = 0; c < dim; c++) xi0[c] = -1.; 78 CHKERRQ(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 79 CHKERRQ(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid)); 80 CHKERRQ(DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **) &vals)); 81 for (c = 0; c < Ncell; ++c) { 82 for (q = 0; q < Nq; ++q) { 83 CHKERRQ(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 CHKERRQ(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 90 CHKERRQ(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid)); 91 CHKERRQ(DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **) &vals)); 92 CHKERRQ(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 CHKERRQ(DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u)); 108 CHKERRQ(VecViewFromOptions(u, NULL, "-f_view")); 109 CHKERRQ(DMGetGlobalVector(dm, &rhs)); 110 CHKERRQ(DMCreateMassMatrix(sw, dm, &mass)); 111 CHKERRQ(MatMult(mass, u, rhs)); 112 CHKERRQ(MatDestroy(&mass)); 113 CHKERRQ(VecDestroy(&u)); 114 115 CHKERRQ(DMGetGlobalVector(dm, &uproj)); 116 CHKERRQ(DMCreateMatrix(dm, &mass)); 117 CHKERRQ(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user)); 118 CHKERRQ(MatViewFromOptions(mass, NULL, "-mass_mat_view")); 119 CHKERRQ(KSPCreate(PETSC_COMM_WORLD, &ksp)); 120 CHKERRQ(KSPSetOperators(ksp, mass, mass)); 121 CHKERRQ(KSPSetFromOptions(ksp)); 122 CHKERRQ(KSPSolve(ksp, rhs, uproj)); 123 CHKERRQ(KSPDestroy(&ksp)); 124 CHKERRQ(PetscObjectSetName((PetscObject) uproj, "Full Projection")); 125 CHKERRQ(VecViewFromOptions(uproj, NULL, "-proj_vec_view")); 126 CHKERRQ(DMComputeL2Diff(dm, 0.0, funcs, NULL, uproj, &error)); 127 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double) error)); 128 CHKERRQ(MatDestroy(&mass)); 129 CHKERRQ(DMRestoreGlobalVector(dm, &rhs)); 130 CHKERRQ(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 PetscErrorCode ierr; 139 140 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 141 comm = PETSC_COMM_WORLD; 142 CHKERRQ(CreateMesh(comm, &dm, &user)); 143 CHKERRQ(CreateParticles(dm, &sw, &user)); 144 CHKERRQ(PetscObjectSetName((PetscObject) dm, "Mesh")); 145 CHKERRQ(DMViewFromOptions(dm, NULL, "-dm_view")); 146 CHKERRQ(DMViewFromOptions(sw, NULL, "-sw_view")); 147 CHKERRQ(TestL2Projection(dm, sw, &user)); 148 CHKERRQ(DMDestroy(&dm)); 149 CHKERRQ(DMDestroy(&sw)); 150 PetscFinalize(); 151 return 0; 152 } 153 154 /*TEST 155 156 test: 157 suffix: proj_0 158 requires: pragmatic 159 TODO: broken 160 args: -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 161 test: 162 suffix: proj_1 163 requires: pragmatic 164 TODO: broken 165 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 166 test: 167 suffix: proj_2 168 requires: pragmatic 169 TODO: broken 170 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 171 test: 172 suffix: proj_3 173 requires: pragmatic 174 TODO: broken 175 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 176 177 TEST*/ 178