1 static char help[] = "Example program demonstrating projection between particle and finite element spaces using OpenMP in 2D cylindrical coordinates\n"; 2 3 #include "petscdmplex.h" 4 #include "petscds.h" 5 #include "petscdmswarm.h" 6 #include "petscksp.h" 7 #include <petsc/private/petscimpl.h> 8 #if defined(PETSC_HAVE_OPENMP) && defined(PETSC_HAVE_THREADSAFETY) 9 #include <omp.h> 10 #endif 11 12 typedef struct { 13 Mat MpTrans; 14 Mat Mp; 15 Vec ff; 16 Vec uu; 17 } MatShellCtx; 18 19 PetscErrorCode MatMultMtM_SeqAIJ(Mat MtM,Vec xx,Vec yy) 20 { 21 MatShellCtx *matshellctx; 22 PetscErrorCode ierr; 23 24 PetscFunctionBeginUser; 25 ierr = MatShellGetContext(MtM,&matshellctx);CHKERRQ(ierr); 26 if (!matshellctx) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "No context"); 27 ierr = MatMult(matshellctx->Mp, xx, matshellctx->ff);CHKERRQ(ierr); 28 ierr = MatMult(matshellctx->MpTrans, matshellctx->ff, yy);CHKERRQ(ierr); 29 PetscFunctionReturn(0); 30 } 31 32 PetscErrorCode MatMultAddMtM_SeqAIJ(Mat MtM,Vec xx, Vec yy, Vec zz) 33 { 34 MatShellCtx *matshellctx; 35 PetscErrorCode ierr; 36 37 PetscFunctionBeginUser; 38 ierr = MatShellGetContext(MtM,&matshellctx);CHKERRQ(ierr); 39 if (!matshellctx) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "No context"); 40 ierr = MatMult(matshellctx->Mp, xx, matshellctx->ff);CHKERRQ(ierr); 41 ierr = MatMultAdd(matshellctx->MpTrans, matshellctx->ff, yy, zz);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 PetscErrorCode createSwarm(const DM dm, DM *sw) 46 { 47 PetscErrorCode ierr; 48 PetscInt Nc = 1, dim = 2; 49 50 PetscFunctionBeginUser; 51 ierr = DMCreate(PETSC_COMM_SELF, sw);CHKERRQ(ierr); 52 ierr = DMSetType(*sw, DMSWARM);CHKERRQ(ierr); 53 ierr = DMSetDimension(*sw, dim);CHKERRQ(ierr); 54 ierr = DMSwarmSetType(*sw, DMSWARM_PIC);CHKERRQ(ierr); 55 ierr = DMSwarmSetCellDM(*sw, dm);CHKERRQ(ierr); 56 ierr = DMSwarmRegisterPetscDatatypeField(*sw, "w_q", Nc, PETSC_SCALAR);CHKERRQ(ierr); 57 ierr = DMSwarmFinalizeFieldRegister(*sw);CHKERRQ(ierr); 58 ierr = DMSetFromOptions(*sw);CHKERRQ(ierr); 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode gridToParticles(const DM dm, DM sw, PetscReal *moments, Vec rhs, Mat M_p) 63 { 64 PetscBool is_lsqr; 65 KSP ksp; 66 Mat PM_p=NULL,MtM,D; 67 Vec ff; 68 PetscErrorCode ierr; 69 PetscInt Np, timestep = 0, bs, N, M, nzl; 70 PetscReal time = 0.0; 71 PetscDataType dtype; 72 MatShellCtx *matshellctx; 73 74 PetscFunctionBeginUser; 75 ierr = KSPCreate(PETSC_COMM_SELF, &ksp);CHKERRQ(ierr); 76 ierr = KSPSetOptionsPrefix(ksp, "ftop_");CHKERRQ(ierr); 77 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 78 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPLSQR,&is_lsqr); 79 if (!is_lsqr) { 80 ierr = MatGetLocalSize(M_p, &M, &N);CHKERRQ(ierr); 81 if (N>M) { 82 PC pc; 83 ierr = PetscInfo2(ksp, " M (%D) < M (%D) -- skip revert to lsqr\n",M,N);CHKERRQ(ierr); 84 is_lsqr = PETSC_TRUE; 85 ierr = KSPSetType(ksp,KSPLSQR);CHKERRQ(ierr); 86 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 87 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); // could put in better solver -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero 88 } else { 89 ierr = PetscNew(&matshellctx);CHKERRQ(ierr); 90 ierr = MatCreateShell(PetscObjectComm((PetscObject)dm),N,N,PETSC_DECIDE,PETSC_DECIDE,matshellctx,&MtM);CHKERRQ(ierr); 91 ierr = MatTranspose(M_p,MAT_INITIAL_MATRIX,&matshellctx->MpTrans);CHKERRQ(ierr); 92 matshellctx->Mp = M_p; 93 ierr = MatShellSetOperation(MtM, MATOP_MULT, (void (*)(void))MatMultMtM_SeqAIJ);CHKERRQ(ierr); 94 ierr = MatShellSetOperation(MtM, MATOP_MULT_ADD, (void (*)(void))MatMultAddMtM_SeqAIJ);CHKERRQ(ierr); 95 ierr = MatCreateVecs(M_p,&matshellctx->uu,&matshellctx->ff);CHKERRQ(ierr); 96 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,1,NULL,&D);CHKERRQ(ierr); 97 for (int i=0 ; i<N ; i++) { 98 const PetscScalar *vals; 99 const PetscInt *cols; 100 PetscScalar dot = 0; 101 ierr = MatGetRow(matshellctx->MpTrans,i,&nzl,&cols,&vals);CHKERRQ(ierr); 102 for (int ii=0 ; ii<nzl ; ii++) dot += PetscSqr(vals[ii]); 103 if (dot==0.0) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Row %D is empty", i); 104 ierr = MatSetValue(D,i,i,dot,INSERT_VALUES); 105 } 106 ierr = MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 107 ierr = MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 108 PetscInfo2(M_p,"createMtMKSP Have %D eqs, nzl = %D\n",N,nzl); 109 ierr = KSPSetOperators(ksp, MtM, D);CHKERRQ(ierr); 110 ierr = MatViewFromOptions(D,NULL,"-ftop2_D_mat_view");CHKERRQ(ierr); 111 ierr = MatViewFromOptions(M_p,NULL,"-ftop2_Mp_mat_view");CHKERRQ(ierr); 112 ierr = MatViewFromOptions(matshellctx->MpTrans,NULL,"-ftop2_MpTranspose_mat_view");CHKERRQ(ierr); 113 } 114 } 115 if (is_lsqr) { 116 PC pc; 117 PetscBool is_bjac; 118 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 119 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&is_bjac); 120 if (is_bjac) { 121 ierr = DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);CHKERRQ(ierr); 122 ierr = KSPSetOperators(ksp, M_p, PM_p);CHKERRQ(ierr); 123 } else { 124 ierr = KSPSetOperators(ksp, M_p, M_p);CHKERRQ(ierr); 125 } 126 } 127 ierr = DMSwarmCreateGlobalVectorFromField(sw, "w_q", &ff);CHKERRQ(ierr); // this grabs access !!!!! 128 if (!is_lsqr) { 129 ierr = KSPSolve(ksp, rhs, matshellctx->uu);CHKERRQ(ierr); 130 ierr = MatMult(M_p, matshellctx->uu, ff);CHKERRQ(ierr); 131 ierr = MatDestroy(&matshellctx->MpTrans);CHKERRQ(ierr); 132 ierr = VecDestroy(&matshellctx->ff);CHKERRQ(ierr); 133 ierr = VecDestroy(&matshellctx->uu);CHKERRQ(ierr); 134 ierr = MatDestroy(&D);CHKERRQ(ierr); 135 ierr = MatDestroy(&MtM);CHKERRQ(ierr); 136 ierr = PetscFree(matshellctx);CHKERRQ(ierr); 137 } else { 138 ierr = KSPSolveTranspose(ksp, rhs, ff);CHKERRQ(ierr); 139 } 140 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 141 /* Visualize particle field */ 142 ierr = DMSetOutputSequenceNumber(sw, timestep, time);CHKERRQ(ierr); 143 ierr = VecViewFromOptions(ff, NULL, "-weights_view");CHKERRQ(ierr); 144 ierr = DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &ff);CHKERRQ(ierr); 145 146 /* compute energy */ 147 if (moments) { 148 PetscReal *wq, *coords; 149 ierr = DMSwarmGetLocalSize(sw,&Np);CHKERRQ(ierr); 150 ierr = DMSwarmGetField(sw, "w_q", &bs, &dtype, (void**)&wq);CHKERRQ(ierr); 151 ierr = DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void**)&coords);CHKERRQ(ierr); 152 moments[0] = moments[1] = moments[2] = 0; 153 for (int p=0;p<Np;p++) { 154 moments[0] += wq[p]; 155 moments[1] += wq[p] * coords[p*2+0]; // x-momentum 156 moments[2] += wq[p] * (PetscSqr(coords[p*2+0])+PetscSqr(coords[p*2+1])); 157 } 158 ierr = DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void**)&coords);CHKERRQ(ierr); 159 ierr = DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void**)&wq);CHKERRQ(ierr); 160 } 161 ierr = MatDestroy(&PM_p); 162 163 PetscFunctionReturn(0); 164 } 165 166 PetscErrorCode particlesToGrid(const DM dm, DM sw, const PetscInt Np, const PetscInt a_tid, const PetscInt dim, const PetscInt target, 167 const PetscReal xx[], const PetscReal yy[], const PetscReal a_wp[], Vec rho, Mat *Mp_out) 168 { 169 170 PetscBool removePoints = PETSC_TRUE; 171 PetscReal *wq, *coords; 172 PetscDataType dtype; 173 Mat M_p; 174 Vec ff; 175 PetscErrorCode ierr; 176 PetscInt bs,p,zero=0; 177 178 PetscFunctionBeginUser; 179 ierr = DMSwarmSetLocalSizes(sw, Np, zero);CHKERRQ(ierr); 180 ierr = DMSwarmGetField(sw, "w_q", &bs, &dtype, (void**)&wq);CHKERRQ(ierr); 181 ierr = DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void**)&coords);CHKERRQ(ierr); 182 for (p=0;p<Np;p++) { 183 coords[p*2+0] = xx[p]; 184 coords[p*2+1] = yy[p]; 185 wq[p] = a_wp[p]; 186 } 187 ierr = DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void**)&coords);CHKERRQ(ierr); 188 ierr = DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void**)&wq);CHKERRQ(ierr); 189 ierr = DMSwarmMigrate(sw, removePoints);CHKERRQ(ierr); 190 ierr = PetscObjectSetName((PetscObject)sw, "Particle Grid");CHKERRQ(ierr); 191 if (a_tid==target) {ierr = DMViewFromOptions(sw, NULL, "-swarm_view");CHKERRQ(ierr);} 192 193 /* Project particles to field */ 194 /* This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE */ 195 ierr = DMCreateMassMatrix(sw, dm, &M_p);CHKERRQ(ierr); 196 ierr = PetscObjectSetName((PetscObject)rho, "rho");CHKERRQ(ierr); 197 198 ierr = DMSwarmCreateGlobalVectorFromField(sw, "w_q", &ff);CHKERRQ(ierr); // this grabs access !!!!! 199 ierr = PetscObjectSetName((PetscObject)ff, "weights");CHKERRQ(ierr); 200 ierr = MatMultTranspose(M_p, ff, rho);CHKERRQ(ierr); 201 ierr = DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &ff);CHKERRQ(ierr); 202 203 /* Visualize mesh field */ 204 if (a_tid==target) {ierr = VecViewFromOptions(rho, NULL, "-rho_view");CHKERRQ(ierr);} 205 // output 206 *Mp_out = M_p; 207 208 PetscFunctionReturn(0); 209 } 210 static PetscErrorCode maxwellian(PetscInt dim, const PetscReal x[], PetscReal kt_m, PetscReal n, PetscScalar *u) 211 { 212 PetscInt i; 213 PetscReal v2 = 0, theta = 2*kt_m; /* theta = 2kT/mc^2 */ 214 215 PetscFunctionBegin; 216 /* compute the exponents, v^2 */ 217 for (i = 0; i < dim; ++i) v2 += x[i]*x[i]; 218 /* evaluate the Maxwellian */ 219 u[0] = n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta)) * 2.*PETSC_PI*x[1]; // radial term for 2D axi-sym. 220 PetscFunctionReturn(0); 221 } 222 #define NUM_SOLVE_LOOPS 100 223 #define MAX_NUM_THRDS 12 224 PetscErrorCode go() 225 { 226 DM dm_t[MAX_NUM_THRDS], sw_t[MAX_NUM_THRDS]; 227 PetscFE fe; 228 PetscInt dim = 2, Nc = 1, timestep = 0, i, faces[3]; 229 PetscInt Np[2] = {10,10}, Np2[2], field = 0, target = 0, Np_t[MAX_NUM_THRDS]; 230 PetscReal time = 0.0, moments_0[3], moments_1[3], vol; 231 PetscReal lo[3] = {-5,0,-5}, hi[3] = {5,5,5}, h[3], hp[3], *xx_t[MAX_NUM_THRDS], *yy_t[MAX_NUM_THRDS], *wp_t[MAX_NUM_THRDS], solve_time = 0; 232 Vec rho_t[MAX_NUM_THRDS], rhs_t[MAX_NUM_THRDS]; 233 Mat M_p_t[MAX_NUM_THRDS]; 234 PetscErrorCode ierr; 235 #if defined PETSC_USE_LOG 236 PetscLogStage stage; 237 PetscLogEvent swarm_create_ev, solve_ev, solve_loop_ev; 238 #endif 239 #if defined(PETSC_HAVE_OPENMP) && defined(PETSC_HAVE_THREADSAFETY) 240 PetscInt numthreads = PetscNumOMPThreads; 241 double starttime, endtime; 242 #else 243 PetscInt numthreads = 1; 244 #endif 245 246 PetscFunctionBeginUser; 247 #if defined(PETSC_HAVE_OPENMP) && defined(PETSC_HAVE_THREADSAFETY) 248 if (numthreads>MAX_NUM_THRDS) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Too many threads %D > %D", numthreads, MAX_NUM_THRDS); 249 if (numthreads<=0) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "No threads %D > %D ", numthreads, MAX_NUM_THRDS); 250 #endif 251 if (target >= numthreads) target = numthreads-1; 252 ierr = PetscLogEventRegister("Create Swarm", DM_CLASSID, &swarm_create_ev);CHKERRQ(ierr); 253 ierr = PetscLogEventRegister("Single solve", DM_CLASSID, &solve_ev);CHKERRQ(ierr); 254 ierr = PetscLogEventRegister("Solve loop", DM_CLASSID, &solve_loop_ev);CHKERRQ(ierr); 255 ierr = PetscLogStageRegister("Solve", &stage);CHKERRQ(ierr); 256 i = dim; 257 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &i, NULL);CHKERRQ(ierr); 258 i = dim; 259 ierr = PetscOptionsGetIntArray(NULL, NULL, "-np", Np, &i, NULL);CHKERRQ(ierr); 260 /* Create thread meshes */ 261 for (int tid=0; tid<numthreads; tid++) { 262 // setup mesh dm_t, could use PETSc's Landau create velocity space mesh here to get dm_t[tid] 263 ierr = DMCreate(PETSC_COMM_SELF, &dm_t[tid]);CHKERRQ(ierr); 264 ierr = DMSetType(dm_t[tid], DMPLEX);CHKERRQ(ierr); 265 ierr = DMSetFromOptions(dm_t[tid]);CHKERRQ(ierr); 266 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, "", PETSC_DECIDE, &fe);CHKERRQ(ierr); 267 ierr = PetscFESetFromOptions(fe);CHKERRQ(ierr); 268 ierr = PetscObjectSetName((PetscObject)fe, "fe");CHKERRQ(ierr); 269 ierr = DMSetField(dm_t[tid], field, NULL, (PetscObject)fe);CHKERRQ(ierr); 270 ierr = DMCreateDS(dm_t[tid]);CHKERRQ(ierr); 271 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 272 // helper vectors 273 ierr = DMSetOutputSequenceNumber(dm_t[tid], timestep, time);CHKERRQ(ierr); // not used 274 ierr = DMCreateGlobalVector(dm_t[tid], &rho_t[tid]);CHKERRQ(ierr); 275 ierr = DMCreateGlobalVector(dm_t[tid], &rhs_t[tid]);CHKERRQ(ierr); 276 // this mimics application code 277 ierr = DMGetBoundingBox(dm_t[tid], lo, hi);CHKERRQ(ierr); 278 if (tid==target) { 279 ierr = DMViewFromOptions(dm_t[tid], NULL, "-dm_view");CHKERRQ(ierr); 280 for (i=0,vol=1;i<dim;i++) { 281 h[i] = (hi[i] - lo[i])/faces[i]; 282 hp[i] = (hi[i] - lo[i])/Np[i]; 283 vol *= (hi[i] - lo[i]); 284 ierr = PetscInfo5(dm_t[tid]," lo = %g hi = %g n = %D h = %g hp = %g\n",lo[i],hi[i],faces[i],h[i],hp[i]);CHKERRQ(ierr); 285 } 286 } 287 } 288 // prepare particle data for problems. This mimics application code 289 ierr = PetscLogEventBegin(swarm_create_ev,0,0,0,0);CHKERRQ(ierr); 290 Np2[0] = Np[0]; Np2[1] = Np[1]; 291 for (int tid=0; tid<numthreads; tid++) { // change size of particle list a little 292 Np_t[tid] = Np2[0]*Np2[1]; 293 ierr = PetscMalloc3(Np_t[tid],&xx_t[tid],Np_t[tid],&yy_t[tid],Np_t[tid],&wp_t[tid]);CHKERRQ(ierr); 294 if (tid==target) {moments_0[0] = moments_0[1] = moments_0[2] = 0;} 295 for (int pi=0, pp=0;pi<Np2[0];pi++) { 296 for (int pj=0;pj<Np2[1];pj++,pp++) { 297 xx_t[tid][pp] = lo[0] + hp[0]/2. + pi*hp[0]; 298 yy_t[tid][pp] = lo[1] + hp[1]/2. + pj*hp[1]; 299 { 300 PetscReal x[] = {xx_t[tid][pp],yy_t[tid][pp]}; 301 ierr = maxwellian(2, x, 1.0, vol/(PetscReal)Np_t[tid], &wp_t[tid][pp]); 302 } 303 if (tid==target) { //energy_0 += wp_t[tid][pp]*(PetscSqr(xx_t[tid][pp])+PetscSqr(yy_t[tid][pp])); 304 moments_0[0] += wp_t[tid][pp]; 305 moments_0[1] += wp_t[tid][pp] * xx_t[tid][pp]; // x-momentum 306 moments_0[2] += wp_t[tid][pp] * (PetscSqr(xx_t[tid][pp]) + PetscSqr(yy_t[tid][pp])); 307 } 308 } 309 } 310 Np2[0]++; Np2[1]++; 311 } 312 ierr = PetscLogEventEnd(swarm_create_ev,0,0,0,0);CHKERRQ(ierr); 313 ierr = PetscLogEventBegin(solve_ev,0,0,0,0);CHKERRQ(ierr); 314 /* Create particle swarm */ 315 PetscPragmaOMP(parallel for) 316 for (int tid=0; tid<numthreads; tid++) { 317 PetscErrorCode ierr_t; 318 ierr_t = createSwarm(dm_t[tid], &sw_t[tid]); 319 if (ierr_t) ierr = ierr_t; 320 } 321 CHKERRQ(ierr); 322 PetscPragmaOMP(parallel for) 323 for (int tid=0; tid<numthreads; tid++) { 324 PetscErrorCode ierr_t; 325 ierr_t = particlesToGrid(dm_t[tid], sw_t[tid], Np_t[tid], tid, dim, target, xx_t[tid], yy_t[tid], wp_t[tid], rho_t[tid], &M_p_t[tid]); 326 if (ierr_t) ierr = ierr_t; 327 } 328 CHKERRQ(ierr); 329 /* Project field to particles */ 330 /* This gives f_p = M_p^+ M f */ 331 PetscPragmaOMP(parallel for) 332 for (int tid=0; tid<numthreads; tid++) { 333 PetscErrorCode ierr_t; 334 ierr_t = VecCopy(rho_t[tid], rhs_t[tid]); /* Identity: M^1 M rho */ 335 if (ierr_t) ierr = ierr_t; 336 } 337 CHKERRQ(ierr); 338 PetscPragmaOMP(parallel for) 339 for (int tid=0; tid<numthreads; tid++) { 340 PetscErrorCode ierr_t; 341 ierr_t = gridToParticles(dm_t[tid], sw_t[tid], (tid==target) ? moments_1 : NULL, rhs_t[tid], M_p_t[tid]); 342 if (ierr_t) ierr = ierr_t; 343 } 344 CHKERRQ(ierr); 345 /* Cleanup */ 346 for (int tid=0; tid<numthreads; tid++) { 347 ierr = MatDestroy(&M_p_t[tid]);CHKERRQ(ierr); 348 ierr = DMDestroy(&sw_t[tid]);CHKERRQ(ierr); 349 } 350 ierr = PetscLogEventEnd(solve_ev,0,0,0,0);CHKERRQ(ierr); 351 /* for timing */ 352 ierr = PetscOptionsClearValue(NULL,"-ftop_ksp_converged_reason");CHKERRQ(ierr); 353 ierr = PetscOptionsClearValue(NULL,"-ftop_ksp_monitor");CHKERRQ(ierr); 354 ierr = PetscOptionsClearValue(NULL,"-ftop_ksp_view");CHKERRQ(ierr); 355 ierr = PetscOptionsClearValue(NULL,"-info");CHKERRQ(ierr); 356 // repeat solve many times to get warmed up data 357 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 358 #if defined(PETSC_HAVE_OPENMP) && defined(PETSC_HAVE_THREADSAFETY) 359 starttime = MPI_Wtime(); 360 #endif 361 ierr = PetscLogEventBegin(solve_loop_ev,0,0,0,0);CHKERRQ(ierr); 362 for (int d=0; d<NUM_SOLVE_LOOPS; d++) { 363 /* Create particle swarm */ 364 PetscPragmaOMP(parallel for) 365 for (int tid=0; tid<numthreads; tid++) { 366 PetscErrorCode ierr_t; 367 ierr_t = createSwarm(dm_t[tid], &sw_t[tid]); 368 if (ierr_t) ierr = ierr_t; 369 } 370 CHKERRQ(ierr); 371 PetscPragmaOMP(parallel for) 372 for (int tid=0; tid<numthreads; tid++) { 373 PetscErrorCode ierr_t; 374 ierr_t = particlesToGrid(dm_t[tid], sw_t[tid], Np_t[tid], tid, dim, target, xx_t[tid], yy_t[tid], wp_t[tid], rho_t[tid], &M_p_t[tid]); 375 if (ierr_t) ierr = ierr_t; 376 } 377 CHKERRQ(ierr); 378 PetscPragmaOMP(parallel for) 379 for (int tid=0; tid<numthreads; tid++) { 380 PetscErrorCode ierr_t; 381 ierr_t = VecCopy(rho_t[tid], rhs_t[tid]); /* Identity: M^1 M rho */ 382 if (ierr_t) ierr = ierr_t; 383 } 384 CHKERRQ(ierr); 385 PetscPragmaOMP(parallel for) 386 for (int tid=0; tid<numthreads; tid++) { 387 PetscErrorCode ierr_t; 388 ierr_t = gridToParticles(dm_t[tid], sw_t[tid], NULL, rhs_t[tid], M_p_t[tid]); 389 if (ierr_t) ierr = ierr_t; 390 } 391 CHKERRQ(ierr); 392 /* Cleanup */ 393 for (int tid=0; tid<numthreads; tid++) { 394 ierr = MatDestroy(&M_p_t[tid]);CHKERRQ(ierr); 395 ierr = DMDestroy(&sw_t[tid]);CHKERRQ(ierr); 396 } 397 } 398 ierr = PetscLogEventEnd(solve_loop_ev,0,0,0,0);CHKERRQ(ierr); 399 #if defined(PETSC_HAVE_OPENMP) && defined(PETSC_HAVE_THREADSAFETY) 400 endtime = MPI_Wtime(); 401 solve_time += (endtime - starttime); 402 #endif 403 ierr = PetscLogStagePop();CHKERRQ(ierr); 404 // 405 ierr = PetscPrintf(PETSC_COMM_SELF,"Total number density: %20.12e (%20.12e); x-momentum = %g (%g); energy = %g error = %e, %D particles. Use %D threads, Solve time: %g\n", moments_1[0], moments_0[0], moments_1[1], moments_0[1], moments_1[2], (moments_1[2]-moments_0[2])/moments_0[2],Np[0]*Np[1],numthreads,solve_time);CHKERRQ(ierr); 406 /* Cleanup */ 407 for (int tid=0; tid<numthreads; tid++) { 408 ierr = VecDestroy(&rho_t[tid]);CHKERRQ(ierr); 409 ierr = VecDestroy(&rhs_t[tid]);CHKERRQ(ierr); 410 ierr = DMDestroy(&dm_t[tid]);CHKERRQ(ierr); 411 ierr = PetscFree3(xx_t[tid],yy_t[tid],wp_t[tid]);CHKERRQ(ierr); 412 } 413 414 PetscFunctionReturn(0); 415 } 416 417 int main(int argc, char **argv) 418 { 419 PetscErrorCode ierr; 420 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 421 ierr = go();CHKERRQ(ierr); 422 ierr = PetscFinalize(); 423 return ierr; 424 } 425 426 /*TEST 427 428 build: 429 requires: !complex 430 431 test: 432 suffix: 0 433 requires: double triangle 434 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -petscspace_degree 2 -ftop_ksp_type lsqr -ftop_pc_type none -dm_view -ftop_ksp_converged_reason -ftop_ksp_rtol 1.e-14 435 filter: grep -v DM_ | grep -v atomic 436 437 test: 438 suffix: 1 439 requires: double triangle 440 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -petscspace_degree 2 -dm_plex_hash_location -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero -dm_view -ftop_ksp_converged_reason -ftop_ksp_rtol 1.e-14 441 filter: grep -v DM_ | grep -v atomic 442 443 test: 444 suffix: 2 445 requires: double triangle 446 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -petscspace_degree 2 -dm_plex_hash_location -ftop_ksp_type cg -ftop_pc_type jacobi -dm_view -ftop_ksp_converged_reason -ftop_ksp_rtol 1.e-14 447 filter: grep -v DM_ | grep -v atomic 448 449 TEST*/ 450