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