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 dim; /* The topological mesh dimension */ 12 PetscBool simplex; /* Flag for simplices or tensor cells */ 13 char mshNam[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */ 14 PetscInt nbrVerEdge; /* Number of vertices per edge if unit square/cube generated */ 15 } AppCtx; 16 17 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 18 { 19 PetscErrorCode ierr; 20 21 PetscFunctionBeginUser; 22 options->dim = 2; 23 options->simplex = PETSC_TRUE; 24 ierr = PetscStrcpy(options->mshNam, "");CHKERRQ(ierr); 25 26 ierr = PetscOptionsBegin(comm, "", "Meshing Adaptation Options", "DMPLEX");CHKERRQ(ierr); 27 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex1.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 28 ierr = PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex1.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 29 ierr = PetscOptionsString("-msh", "Name of the mesh filename if any", "ex1.c", options->mshNam, options->mshNam, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 30 ierr = PetscOptionsInt("-nbrVerEdge", "Number of vertices per edge if unit square/cube generated", "ex1.c", options->nbrVerEdge, &options->nbrVerEdge, NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsEnd(); 32 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 37 { 38 PetscBool flag; 39 PetscErrorCode ierr; 40 41 PetscFunctionBeginUser; 42 ierr = PetscStrcmp(user->mshNam, "", &flag);CHKERRQ(ierr); 43 if (flag) { 44 PetscInt faces[3]; 45 46 faces[0] = user->nbrVerEdge-1; faces[1] = user->nbrVerEdge-1; faces[2] = user->nbrVerEdge-1; 47 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, faces, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 48 } else { 49 ierr = DMPlexCreateFromFile(comm, user->mshNam, PETSC_TRUE, dm);CHKERRQ(ierr); 50 ierr = DMGetDimension(*dm, &user->dim);CHKERRQ(ierr); 51 } 52 { 53 DM distributedMesh = NULL; 54 55 /* Distribute mesh over processes */ 56 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 57 if (distributedMesh) { 58 ierr = DMDestroy(dm);CHKERRQ(ierr); 59 *dm = distributedMesh; 60 } 61 } 62 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 67 { 68 PetscInt d; 69 u[0] = 0.0; 70 for (d = 0; d < dim; ++d) u[0] += x[d]; 71 return 0; 72 } 73 74 static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, 75 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 76 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 77 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 78 { 79 g0[0] = 1.0; 80 } 81 82 static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user) 83 { 84 PetscDS prob; 85 PetscFE fe; 86 PetscQuadrature quad; 87 PetscScalar *vals; 88 PetscReal *v0, *J, *invJ, detJ, *coords, *xi0; 89 PetscInt *cellid; 90 const PetscReal *qpoints; 91 PetscInt Ncell, c, Nq, q, dim; 92 PetscErrorCode ierr; 93 94 PetscFunctionBeginUser; 95 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 96 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, user->simplex, NULL, -1, &fe);CHKERRQ(ierr); 97 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 98 ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe);CHKERRQ(ierr); 99 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 100 ierr = PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL);CHKERRQ(ierr); 101 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);CHKERRQ(ierr); 102 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 103 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 104 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);CHKERRQ(ierr); 105 106 ierr = DMCreate(PetscObjectComm((PetscObject) dm), sw);CHKERRQ(ierr); 107 ierr = DMSetType(*sw, DMSWARM);CHKERRQ(ierr); 108 ierr = DMSetDimension(*sw, dim);CHKERRQ(ierr); 109 110 ierr = DMSwarmSetType(*sw, DMSWARM_PIC);CHKERRQ(ierr); 111 ierr = DMSwarmSetCellDM(*sw, dm);CHKERRQ(ierr); 112 ierr = DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_SCALAR);CHKERRQ(ierr); 113 ierr = DMSwarmFinalizeFieldRegister(*sw);CHKERRQ(ierr); 114 ierr = DMSwarmSetLocalSizes(*sw, Ncell * Nq, 0);CHKERRQ(ierr); 115 ierr = DMSetFromOptions(*sw);CHKERRQ(ierr); 116 117 ierr = PetscMalloc4(dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);CHKERRQ(ierr); 118 for (c = 0; c < dim; c++) xi0[c] = -1.; 119 ierr = DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 120 ierr = DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 121 ierr = DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **) &vals);CHKERRQ(ierr); 122 for (c = 0; c < Ncell; ++c) { 123 for (q = 0; q < Nq; ++q) { 124 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 125 cellid[c*Nq + q] = c; 126 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], &coords[(c*Nq + q)*dim]); 127 linear(dim, 0.0, &coords[(c*Nq + q)*dim], 1, &vals[c*Nq + q], NULL); 128 } 129 } 130 ierr = DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 131 ierr = DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 132 ierr = DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **) &vals);CHKERRQ(ierr); 133 ierr = PetscFree4(xi0, v0, J, invJ);CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user) 138 { 139 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *); 140 KSP ksp; 141 Mat mass; 142 Vec u, rhs, uproj; 143 PetscReal error; 144 PetscErrorCode ierr; 145 146 PetscFunctionBeginUser; 147 funcs[0] = linear; 148 149 ierr = DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u);CHKERRQ(ierr); 150 ierr = VecViewFromOptions(u, NULL, "-f_view");CHKERRQ(ierr); 151 ierr = DMGetGlobalVector(dm, &rhs);CHKERRQ(ierr); 152 ierr = DMCreateMassMatrix(sw, dm, &mass);CHKERRQ(ierr); 153 ierr = MatMult(mass, u, rhs);CHKERRQ(ierr); 154 ierr = MatDestroy(&mass);CHKERRQ(ierr); 155 ierr = VecDestroy(&u);CHKERRQ(ierr); 156 157 ierr = DMGetGlobalVector(dm, &uproj);CHKERRQ(ierr); 158 ierr = DMCreateMatrix(dm, &mass);CHKERRQ(ierr); 159 ierr = DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user);CHKERRQ(ierr); 160 ierr = MatViewFromOptions(mass, NULL, "-mass_mat_view");CHKERRQ(ierr); 161 ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);CHKERRQ(ierr); 162 ierr = KSPSetOperators(ksp, mass, mass);CHKERRQ(ierr); 163 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 164 ierr = KSPSolve(ksp, rhs, uproj);CHKERRQ(ierr); 165 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 166 ierr = PetscObjectSetName((PetscObject) uproj, "Full Projection");CHKERRQ(ierr); 167 ierr = VecViewFromOptions(uproj, NULL, "-proj_vec_view");CHKERRQ(ierr); 168 ierr = DMComputeL2Diff(dm, 0.0, funcs, NULL, uproj, &error);CHKERRQ(ierr); 169 ierr = PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double) error);CHKERRQ(ierr); 170 ierr = MatDestroy(&mass);CHKERRQ(ierr); 171 ierr = DMRestoreGlobalVector(dm, &rhs);CHKERRQ(ierr); 172 ierr = DMRestoreGlobalVector(dm, &uproj);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 int main (int argc, char * argv[]) { 177 MPI_Comm comm; 178 DM dm, sw; 179 AppCtx user; 180 PetscErrorCode ierr; 181 182 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 183 comm = PETSC_COMM_WORLD; 184 ierr = ProcessOptions(comm, &user);CHKERRQ(ierr); 185 ierr = CreateMesh(comm, &dm, &user);CHKERRQ(ierr); 186 ierr = CreateParticles(dm, &sw, &user);CHKERRQ(ierr); 187 ierr = PetscObjectSetName((PetscObject) dm, "Mesh");CHKERRQ(ierr); 188 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); 189 ierr = DMViewFromOptions(sw, NULL, "-sw_view");CHKERRQ(ierr); 190 ierr = TestL2Projection(dm, sw, &user);CHKERRQ(ierr); 191 ierr = DMDestroy(&dm);CHKERRQ(ierr); 192 ierr = DMDestroy(&sw);CHKERRQ(ierr); 193 PetscFinalize(); 194 return 0; 195 } 196 197 /*TEST 198 199 test: 200 suffix: proj_0 201 requires: pragmatic 202 TODO: broken 203 args: -dim 2 -nbrVerEdge 3 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu 204 test: 205 suffix: proj_1 206 requires: pragmatic 207 TODO: broken 208 args: -dim 2 -simplex 0 -nbrVerEdge 3 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu 209 test: 210 suffix: proj_2 211 requires: pragmatic 212 TODO: broken 213 args: -dim 3 -nbrVerEdge 3 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu 214 test: 215 suffix: proj_3 216 requires: pragmatic 217 TODO: broken 218 args: -dim 2 -simplex 0 -nbrVerEdge 3 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu 219 220 TEST*/ 221