Actual source code: ex30.c
1: static char help[] = "Grid based Landau collision operator with PIC interface with OpenMP setup. (one species per grid)\n";
3: /*
4: Support 2.5V with axisymmetric coordinates
5: - r,z coordinates
6: - Domain and species data input by Landau operator
7: - "radius" for each grid, normalized with electron thermal velocity
8: - Domain: (0,radius) x (-radius,radius), thus first coordinate x[0] is perpendicular velocity and 2pi*x[0] term is added for axisymmetric
9: Supports full 3V
11: */
13: #include <petscdmplex.h>
14: #include <petscds.h>
15: #include <petscdmswarm.h>
16: #include <petscksp.h>
17: #include <petsc/private/petscimpl.h>
18: #if defined(PETSC_HAVE_OPENMP) && defined(PETSC_HAVE_THREADSAFETY)
19: #include <omp.h>
20: #endif
21: #include <petsclandau.h>
22: #include <petscdmcomposite.h>
24: typedef struct {
25: Mat MpTrans;
26: Mat Mp;
27: Vec ff;
28: Vec uu;
29: } MatShellCtx;
31: typedef struct {
32: PetscInt v_target;
33: PetscInt g_target;
34: PetscInt global_vertex_id_0;
35: DM *globSwarmArray;
36: LandauCtx *ctx;
37: DM *grid_dm;
38: Mat *g_Mass;
39: Mat *globMpArray;
40: Vec *globXArray;
41: PetscBool print;
42: PetscBool print_entropy;
43: } PrintCtx;
45: PetscErrorCode MatMultMtM_SeqAIJ(Mat MtM, Vec xx, Vec yy)
46: {
47: MatShellCtx *matshellctx;
49: PetscFunctionBeginUser;
50: PetscCall(MatShellGetContext(MtM, &matshellctx));
51: PetscCheck(matshellctx, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "No context");
52: PetscCall(MatMult(matshellctx->Mp, xx, matshellctx->ff));
53: PetscCall(MatMult(matshellctx->MpTrans, matshellctx->ff, yy));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: PetscErrorCode MatMultAddMtM_SeqAIJ(Mat MtM, Vec xx, Vec yy, Vec zz)
58: {
59: MatShellCtx *matshellctx;
61: PetscFunctionBeginUser;
62: PetscCall(MatShellGetContext(MtM, &matshellctx));
63: PetscCheck(matshellctx, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "No context");
64: PetscCall(MatMult(matshellctx->Mp, xx, matshellctx->ff));
65: PetscCall(MatMultAdd(matshellctx->MpTrans, matshellctx->ff, yy, zz));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: PetscErrorCode createSwarm(const DM dm, PetscInt dim, DM *sw)
70: {
71: PetscInt Nc = 1;
73: PetscFunctionBeginUser;
74: PetscCall(DMCreate(PETSC_COMM_SELF, sw));
75: PetscCall(DMSetType(*sw, DMSWARM));
76: PetscCall(DMSetDimension(*sw, dim));
77: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
78: PetscCall(DMSwarmSetCellDM(*sw, dm));
79: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", Nc, PETSC_REAL));
80: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
81: PetscCall(DMSetFromOptions(*sw));
82: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particle Grid"));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode makeSwarm(DM sw, const PetscInt dim, const PetscInt Np, const PetscReal xx[], const PetscReal yy[], const PetscReal zz[])
87: {
88: PetscReal *coords;
89: PetscDataType dtype;
90: PetscInt bs, p, zero = 0;
92: PetscFunctionBeginUser;
93: PetscCall(DMSwarmSetLocalSizes(sw, Np, zero));
94: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
95: for (p = 0; p < Np; p++) {
96: coords[p * dim + 0] = xx[p];
97: coords[p * dim + 1] = yy[p];
98: if (dim == 3) coords[p * dim + 2] = zz[p];
99: }
100: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
101: PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode createMp(const DM dm, DM sw, Mat *Mp_out)
106: {
107: PetscBool removePoints = PETSC_TRUE;
108: Mat M_p;
110: PetscFunctionBeginUser;
111: // migrate after coords are set
112: PetscCall(DMSwarmMigrate(sw, removePoints));
113: //
114: PetscCall(PetscObjectSetName((PetscObject)sw, "Particle Grid"));
116: /* PetscInt N,*count,nmin=10000,nmax=0,ntot=0; */
117: /* // count */
118: /* PetscCall(DMSwarmCreatePointPerCellCount(sw, &N, &count)); */
119: /* for (int i=0, n; i< N ; i++) { */
120: /* if ((n=count[i]) > nmax) nmax = n; */
121: /* if (n < nmin) nmin = n; */
122: /* PetscCall(PetscInfo(dm, " %d) %d particles\n", i, n)); */
123: /* ntot += n; */
124: /* } */
125: /* PetscCall(PetscFree(count)); */
126: /* PetscCall(PetscInfo(dm, " %" PetscInt_FMT " max particle / cell, and %" PetscInt_FMT " min, ratio = %g, %" PetscInt_FMT " total\n", nmax, nmin, (double)nmax/(double)nmin,ntot)); */
128: /* This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE */
129: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
130: PetscCall(DMViewFromOptions(sw, NULL, "-ex30_sw_view"));
131: // output
132: *Mp_out = M_p;
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: static PetscErrorCode particlesToGrid(const DM dm, DM sw, const PetscInt a_tid, const PetscInt dim, const PetscReal a_wp[], Vec rho, Mat M_p)
137: {
138: PetscReal *wq;
139: PetscDataType dtype;
140: Vec ff;
141: PetscInt bs, p, Np;
143: PetscFunctionBeginUser;
144: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq));
145: PetscCall(DMSwarmGetLocalSize(sw, &Np));
146: for (p = 0; p < Np; p++) wq[p] = a_wp[p];
147: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq));
148: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
149: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &ff));
150: PetscCall(PetscObjectSetName((PetscObject)ff, "weights"));
151: PetscCall(MatMultTranspose(M_p, ff, rho));
152: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &ff));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: //
157: // add grid to arg 'sw.w_q'
158: //
159: PetscErrorCode gridToParticles(const DM dm, DM sw, const Vec rhs, Vec work_ferhs, Mat M_p, Mat Mass)
160: {
161: PetscBool is_lsqr;
162: KSP ksp;
163: Mat PM_p = NULL, MtM, D = NULL;
164: Vec ff;
165: PetscInt N, M, nzl;
166: MatShellCtx *matshellctx = NULL;
167: PC pc;
169: PetscFunctionBeginUser;
170: // 1) apply M in, for Moore-Penrose with mass: Mp (Mp' Mp)^-1 M
171: PetscCall(MatMult(Mass, rhs, work_ferhs));
172: // 2) pseudo-inverse, first part: (Mp' Mp)^-1
173: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
174: PetscCall(KSPSetType(ksp, KSPCG));
175: PetscCall(KSPGetPC(ksp, &pc));
176: PetscCall(PCSetType(pc, PCJACOBI));
177: PetscCall(KSPSetOptionsPrefix(ksp, "ftop_"));
178: PetscCall(KSPSetFromOptions(ksp));
179: PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPLSQR, &is_lsqr));
180: if (!is_lsqr) {
181: PetscCall(MatGetLocalSize(M_p, &M, &N));
182: if (N > M) {
183: PetscCall(PetscInfo(ksp, " M (%" PetscInt_FMT ") < M (%" PetscInt_FMT ") more vertices than particles: revert to lsqr\n", M, N));
184: is_lsqr = PETSC_TRUE;
185: PetscCall(KSPSetType(ksp, KSPLSQR));
186: PetscCall(PCSetType(pc, PCNONE)); // should not happen, but could solve stable (Mp^T Mp), move projection Mp before solve
187: } else {
188: PetscCall(PetscNew(&matshellctx));
189: PetscCall(MatCreateVecs(M_p, &matshellctx->uu, &matshellctx->ff));
190: if (0) {
191: PetscCall(MatTransposeMatMult(M_p, M_p, MAT_INITIAL_MATRIX, 4, &MtM));
192: PetscCall(KSPSetOperators(ksp, MtM, MtM));
193: PetscCall(PetscInfo(M_p, "createMtM KSP with explicit Mp'Mp\n"));
194: PetscCall(MatViewFromOptions(MtM, NULL, "-ftop2_MtM_mat_view"));
195: } else {
196: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)dm), N, N, PETSC_DECIDE, PETSC_DECIDE, matshellctx, &MtM));
197: PetscCall(MatTranspose(M_p, MAT_INITIAL_MATRIX, &matshellctx->MpTrans));
198: matshellctx->Mp = M_p;
199: PetscCall(MatShellSetOperation(MtM, MATOP_MULT, (void (*)(void))MatMultMtM_SeqAIJ));
200: PetscCall(MatShellSetOperation(MtM, MATOP_MULT_ADD, (void (*)(void))MatMultAddMtM_SeqAIJ));
201: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 1, NULL, &D));
202: PetscCall(MatViewFromOptions(matshellctx->MpTrans, NULL, "-ftop2_MpT_mat_view"));
203: for (PetscInt i = 0; i < N; i++) {
204: const PetscScalar *vals;
205: const PetscInt *cols;
206: PetscScalar dot = 0;
207: PetscCall(MatGetRow(matshellctx->MpTrans, i, &nzl, &cols, &vals));
208: for (PetscInt ii = 0; ii < nzl; ii++) dot += PetscSqr(vals[ii]);
209: if (dot < PETSC_MACHINE_EPSILON) {
210: PetscCall(PetscInfo(ksp, "empty row in pseudo-inverse %d\n", (int)i));
211: is_lsqr = PETSC_TRUE; // empty rows
212: PetscCall(KSPSetType(ksp, KSPLSQR));
213: PetscCall(PCSetType(pc, PCNONE)); // should not happen, but could solve stable (Mp Mp^T), move projection Mp before solve
214: // clean up
215: PetscCall(MatDestroy(&matshellctx->MpTrans));
216: PetscCall(VecDestroy(&matshellctx->ff));
217: PetscCall(VecDestroy(&matshellctx->uu));
218: PetscCall(MatDestroy(&D));
219: PetscCall(MatDestroy(&MtM));
220: PetscCall(PetscFree(matshellctx));
221: D = NULL;
222: break;
223: }
224: PetscCall(MatSetValue(D, i, i, dot, INSERT_VALUES));
225: }
226: if (D) {
227: PetscCall(MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY));
228: PetscCall(MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY));
229: PetscCall(PetscInfo(M_p, "createMtMKSP Have %" PetscInt_FMT " eqs, nzl = %" PetscInt_FMT "\n", N, nzl));
230: PetscCall(KSPSetOperators(ksp, MtM, D));
231: PetscCall(MatViewFromOptions(D, NULL, "-ftop2_D_mat_view"));
232: PetscCall(MatViewFromOptions(M_p, NULL, "-ftop2_Mp_mat_view"));
233: PetscCall(MatViewFromOptions(matshellctx->MpTrans, NULL, "-ftop2_MpTranspose_mat_view"));
234: PetscCall(MatViewFromOptions(MtM, NULL, "-ftop2_MtM_mat_view"));
235: }
236: }
237: }
238: }
239: if (is_lsqr) {
240: PC pc2;
241: PetscBool is_bjac;
242: PetscCall(KSPGetPC(ksp, &pc2));
243: PetscCall(PetscObjectTypeCompare((PetscObject)pc2, PCBJACOBI, &is_bjac));
244: if (is_bjac) {
245: PetscCall(DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p));
246: PetscCall(KSPSetOperators(ksp, M_p, PM_p));
247: } else {
248: PetscCall(KSPSetOperators(ksp, M_p, M_p));
249: }
250: PetscCall(MatViewFromOptions(M_p, NULL, "-ftop2_Mp_mat_view"));
251: }
252: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &ff)); // this grabs access
253: if (!is_lsqr) {
254: PetscErrorCode ierr;
255: ierr = KSPSolve(ksp, work_ferhs, matshellctx->uu);
256: if (!ierr) {
257: // 3) with Moore-Penrose apply Mp: M_p (Mp' Mp)^-1 M
258: PetscCall(MatMult(M_p, matshellctx->uu, ff));
259: } else { // failed
260: PC pc2;
261: PetscBool is_bjac;
262: PetscCall(PetscInfo(ksp, "Solver failed, probably singular, try lsqr\n"));
263: PetscCall(KSPReset(ksp));
264: PetscCall(KSPSetType(ksp, KSPLSQR));
265: PetscCall(KSPGetPC(ksp, &pc2));
266: PetscCall(PCSetType(pc2, PCNONE)); // should not happen, but could solve stable (Mp Mp^T), move projection Mp before solve
267: PetscCall(KSPSetOptionsPrefix(ksp, "ftop_"));
268: PetscCall(KSPSetFromOptions(ksp));
269: PetscCall(PetscObjectTypeCompare((PetscObject)pc2, PCBJACOBI, &is_bjac));
270: if (is_bjac) {
271: PetscCall(DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p));
272: PetscCall(KSPSetOperators(ksp, M_p, PM_p));
273: } else {
274: PetscCall(KSPSetOperators(ksp, M_p, M_p));
275: }
276: ierr = KSPSolveTranspose(ksp, work_ferhs, ff);
277: if (ierr) { PetscCheck(!ierr, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "backup LSQR solver failed - need to add N_v > N_p Moore-Penrose pseudo-inverse"); }
278: }
279: if (D) PetscCall(MatDestroy(&D));
280: PetscCall(MatDestroy(&MtM));
281: if (matshellctx->MpTrans) PetscCall(MatDestroy(&matshellctx->MpTrans));
282: PetscCall(VecDestroy(&matshellctx->ff));
283: PetscCall(VecDestroy(&matshellctx->uu));
284: PetscCall(PetscFree(matshellctx));
285: } else {
286: PetscErrorCode ierr;
287: // finally with LSQR apply M_p^\dagger
288: ierr = KSPSolveTranspose(ksp, work_ferhs, ff);
289: if (ierr) { PetscCheck(!ierr, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "backup LSQR solver failed - need to add N_v > N_p Moore-Penrose pseudo-inverse"); }
290: }
291: PetscCall(KSPDestroy(&ksp));
292: PetscCall(MatDestroy(&PM_p));
293: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &ff));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: #define EX30_MAX_NUM_THRDS 12
298: #define EX30_MAX_BATCH_SZ 1024
299: //
300: // add grid to arg 'globSwarmArray[].w_q'
301: //
302: PetscErrorCode gridToParticles_private(DM grid_dm[], DM globSwarmArray[], const PetscInt dim, const PetscInt v_target, const PetscInt numthreads, const PetscInt num_vertices, const PetscInt global_vertex_id, Mat globMpArray[], Mat g_Mass[], Vec t_fhat[][EX30_MAX_NUM_THRDS], PetscReal moments[], Vec globXArray[], LandauCtx *ctx)
303: {
304: PetscErrorCode ierr = (PetscErrorCode)0; // used for inside thread loops
306: PetscFunctionBeginUser;
307: // map back to particles
308: for (PetscInt v_id_0 = 0; v_id_0 < ctx->batch_sz; v_id_0 += numthreads) {
309: PetscCall(PetscInfo(grid_dm[0], "g2p: global batch %" PetscInt_FMT " of %" PetscInt_FMT ", Landau batch %" PetscInt_FMT " of %" PetscInt_FMT ": map back to particles\n", global_vertex_id + 1, num_vertices, v_id_0 + 1, ctx->batch_sz));
310: //PetscPragmaOMP(parallel for)
311: for (PetscInt tid = 0; tid < numthreads; tid++) {
312: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id + v_id;
313: if (glb_v_id < num_vertices) {
314: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
315: PetscErrorCode ierr_t;
316: ierr_t = PetscInfo(grid_dm[0], "gridToParticles: global batch %" PetscInt_FMT ", local batch b=%" PetscInt_FMT ", grid g=%" PetscInt_FMT ", index(b,g) %" PetscInt_FMT "\n", global_vertex_id, v_id, grid, LAND_PACK_IDX(v_id, grid));
317: ierr_t = gridToParticles(grid_dm[grid], globSwarmArray[LAND_PACK_IDX(v_id, grid)], globXArray[LAND_PACK_IDX(v_id, grid)], t_fhat[grid][tid], globMpArray[LAND_PACK_IDX(v_id, grid)], g_Mass[grid]);
318: if (ierr_t) ierr = ierr_t;
319: }
320: }
321: }
322: PetscCheck(!ierr, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in OMP loop. ierr = %d", (int)ierr);
323: /* Get moments */
324: PetscCall(PetscInfo(grid_dm[0], "Cleanup batches %" PetscInt_FMT " to %" PetscInt_FMT "\n", v_id_0, v_id_0 + numthreads));
325: for (PetscInt tid = 0; tid < numthreads; tid++) {
326: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id + v_id;
327: if (glb_v_id == v_target) {
328: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
329: PetscDataType dtype;
330: PetscReal *wp, *coords;
331: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
332: PetscInt npoints, bs = 1;
333: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wp)); // take data out here
334: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
335: PetscCall(DMSwarmGetLocalSize(sw, &npoints));
336: for (PetscInt p = 0; p < npoints; p++) {
337: PetscReal v2 = 0, fact = (dim == 2) ? 2.0 * PETSC_PI * coords[p * dim + 0] : 1, w = fact * wp[p] * ctx->n_0 * ctx->masses[ctx->species_offset[grid]];
338: for (PetscInt i = 0; i < dim; ++i) v2 += PetscSqr(coords[p * dim + i]);
339: moments[0] += w;
340: moments[1] += w * ctx->v_0 * coords[p * dim + 1]; // z-momentum
341: moments[2] += w * 0.5 * ctx->v_0 * ctx->v_0 * v2;
342: }
343: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wp));
344: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
345: }
346: const PetscReal N_inv = 1 / moments[0];
347: PetscCall(PetscInfo(grid_dm[0], "gridToParticles_private [%" PetscInt_FMT "], n = %g\n", v_id, (double)moments[0]));
348: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
349: PetscDataType dtype;
350: PetscReal *wp, *coords;
351: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
352: PetscInt npoints, bs = 1;
353: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wp)); // take data out here
354: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
355: PetscCall(DMSwarmGetLocalSize(sw, &npoints));
356: for (PetscInt p = 0; p < npoints; p++) {
357: const PetscReal fact = dim == 2 ? 2.0 * PETSC_PI * coords[p * dim + 0] : 1, w = fact * wp[p] * ctx->n_0 * ctx->masses[ctx->species_offset[grid]], ww = w * N_inv;
358: if (w > PETSC_REAL_MIN) {
359: moments[3] -= ww * PetscLogReal(ww);
360: PetscCheck(ww < 1 - PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ww (%g) > 1", (double)ww);
361: } else moments[4] -= w; // keep track of density that is lost
362: }
363: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wp));
364: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
365: }
366: }
367: } // thread batch
368: } // batch
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: static void maxwellian(PetscInt dim, const PetscReal x[], PetscReal kt_m, PetscReal n, PetscReal shift, PetscScalar *u)
373: {
374: PetscInt i;
375: PetscReal v2 = 0, theta = 2.0 * kt_m; /* theta = 2kT/mc^2 */
377: if (shift != 0.) {
378: v2 = 0;
379: for (i = 0; i < dim - 1; ++i) v2 += x[i] * x[i];
380: v2 += (x[dim - 1] - shift) * (x[dim - 1] - shift);
381: /* evaluate the shifted Maxwellian */
382: u[0] += n * PetscPowReal(PETSC_PI * theta, -1.5) * (PetscExpReal(-v2 / theta));
383: } else {
384: /* compute the exponents, v^2 */
385: for (i = 0; i < dim; ++i) v2 += x[i] * x[i];
386: /* evaluate the Maxwellian */
387: u[0] += n * PetscPowReal(PETSC_PI * theta, -1.5) * (PetscExpReal(-v2 / theta));
388: }
389: }
391: static PetscErrorCode PostStep(TS ts)
392: {
393: PetscInt n, dim, nDMs, v_id;
394: PetscReal t;
395: LandauCtx *ctx;
396: Vec X;
397: PrintCtx *printCtx;
398: DM pack;
399: PetscReal moments[5], e_grid[LANDAU_MAX_GRIDS];
401: PetscFunctionBeginUser;
402: PetscCall(TSGetApplicationContext(ts, &printCtx));
403: if (!printCtx->print && !printCtx->print_entropy) PetscFunctionReturn(PETSC_SUCCESS);
404: ctx = printCtx->ctx;
405: if (printCtx->v_target < printCtx->global_vertex_id_0 || printCtx->v_target >= printCtx->global_vertex_id_0 + ctx->batch_sz) PetscFunctionReturn(PETSC_SUCCESS);
406: for (PetscInt i = 0; i < 5; i++) moments[i] = 0;
407: for (PetscInt i = 0; i < LANDAU_MAX_GRIDS; i++) e_grid[i] = 0;
408: v_id = printCtx->v_target % ctx->batch_sz;
409: PetscCall(TSGetDM(ts, &pack));
410: PetscCall(DMGetDimension(pack, &dim));
411: PetscCall(DMCompositeGetNumberDM(pack, &nDMs)); // number of vertices * number of grids
412: PetscCall(TSGetSolution(ts, &X));
413: PetscCall(TSGetStepNumber(ts, &n));
414: PetscCall(TSGetTime(ts, &t));
415: PetscCall(DMCompositeGetAccessArray(pack, X, nDMs, NULL, printCtx->globXArray));
416: if (printCtx->print_entropy && printCtx->v_target >= 0 && 0) {
417: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
418: PetscDataType dtype;
419: PetscReal *wp, *coords;
420: DM sw = printCtx->globSwarmArray[LAND_PACK_IDX(v_id, grid)];
421: Vec work, subX = printCtx->globXArray[LAND_PACK_IDX(v_id, grid)];
422: PetscInt bs, NN;
423: // C-G moments
424: PetscCall(VecDuplicate(subX, &work));
425: PetscCall(gridToParticles(printCtx->grid_dm[grid], sw, subX, work, printCtx->globMpArray[LAND_PACK_IDX(v_id, grid)], printCtx->g_Mass[grid]));
426: PetscCall(VecDestroy(&work));
427: // moments
428: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
429: PetscCall(DMSwarmGetLocalSize(sw, &NN));
430: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wp));
431: for (PetscInt pp = 0; pp < NN; pp++) {
432: PetscReal v2 = 0, fact = (dim == 2) ? 2.0 * PETSC_PI * coords[pp * dim + 0] : 1, w = fact * wp[pp] * ctx->n_0 * ctx->masses[ctx->species_offset[grid]];
433: for (PetscInt i = 0; i < dim; ++i) v2 += PetscSqr(coords[pp * dim + i]);
434: moments[0] += w;
435: moments[1] += w * ctx->v_0 * coords[pp * dim + 1]; // z-momentum
436: moments[2] += w * 0.5 * ctx->v_0 * ctx->v_0 * v2;
437: e_grid[grid] += w * 0.5 * ctx->v_0 * ctx->v_0 * v2;
438: }
439: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
440: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wp));
441: }
442: // entropy
443: const PetscReal N_inv = 1 / moments[0];
444: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
445: PetscDataType dtype;
446: PetscReal *wp, *coords;
447: DM sw = printCtx->globSwarmArray[LAND_PACK_IDX(v_id, grid)];
448: PetscInt bs, NN;
449: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
450: PetscCall(DMSwarmGetLocalSize(sw, &NN));
451: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wp));
452: for (PetscInt pp = 0; pp < NN; pp++) {
453: PetscReal fact = (dim == 2) ? 2.0 * PETSC_PI * coords[pp * dim + 0] : 1, w = fact * wp[pp] * ctx->n_0 * ctx->masses[ctx->species_offset[grid]], ww = w * N_inv;
454: if (w > PETSC_REAL_MIN) {
455: moments[3] -= ww * PetscLogReal(ww);
456: } else moments[4] -= w;
457: }
458: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
459: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wp));
460: }
461: PetscCall(PetscInfo(X, "%4d) time %e, Landau particle moments: 0: %18.12e 1: %19.12e 2: %18.12e entropy: %e loss %e. energy = %e + %e + %e\n", (int)n, (double)t, (double)moments[0], (double)moments[1], (double)moments[2], (double)moments[3], (double)(moments[4] / moments[0]), (double)e_grid[0], (double)e_grid[1], (double)e_grid[2]));
462: }
463: if (printCtx->print && printCtx->g_target >= 0) {
464: PetscInt grid = printCtx->g_target, id;
465: static PetscReal last_t = -100000, period = .5;
466: if (last_t == -100000) last_t = -period + t;
467: if (t >= last_t + period) {
468: last_t = t;
469: PetscCall(DMGetOutputSequenceNumber(ctx->plex[grid], &id, NULL));
470: PetscCall(DMSetOutputSequenceNumber(ctx->plex[grid], id + 1, t));
471: PetscCall(VecViewFromOptions(printCtx->globXArray[LAND_PACK_IDX(v_id % ctx->batch_sz, grid)], NULL, "-ex30_vec_view"));
472: if (ctx->num_grids > grid + 1) {
473: PetscCall(DMSetOutputSequenceNumber(ctx->plex[grid + 1], id + 1, t));
474: PetscCall(VecViewFromOptions(printCtx->globXArray[LAND_PACK_IDX(v_id % ctx->batch_sz, grid + 1)], NULL, "-ex30_vec_view2"));
475: }
476: PetscCall(PetscInfo(X, "%4d) time %e View\n", (int)n, (double)t));
477: }
478: }
479: PetscCall(DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, printCtx->globXArray));
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: PetscErrorCode go(TS ts, Vec X, const PetscInt num_vertices, const PetscInt a_Np, const PetscInt dim, const PetscInt v_target, const PetscInt g_target, PetscReal shift, PetscBool use_uniform_particle_grid)
484: {
485: DM pack, *globSwarmArray, grid_dm[LANDAU_MAX_GRIDS];
486: Mat *globMpArray, g_Mass[LANDAU_MAX_GRIDS];
487: KSP t_ksp[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS];
488: Vec t_fhat[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS];
489: PetscInt nDMs;
490: PetscErrorCode ierr = (PetscErrorCode)0; // used for inside thread loops
491: #if defined(PETSC_HAVE_OPENMP) && defined(PETSC_HAVE_THREADSAFETY)
492: PetscInt numthreads = PetscNumOMPThreads;
493: #else
494: PetscInt numthreads = 1;
495: #endif
496: LandauCtx *ctx;
497: Vec *globXArray;
498: PetscReal moments_0[5], moments_1a[5], moments_1b[5], dt_init;
499: PrintCtx *printCtx;
501: PetscFunctionBeginUser;
502: PetscCheck(numthreads <= EX30_MAX_NUM_THRDS, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Too many threads %" PetscInt_FMT " > %d", numthreads, EX30_MAX_NUM_THRDS);
503: PetscCheck(numthreads > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number threads %" PetscInt_FMT " > %d", numthreads, EX30_MAX_NUM_THRDS);
504: PetscCall(TSGetDM(ts, &pack));
505: PetscCall(DMGetApplicationContext(pack, &ctx));
506: PetscCheck(ctx->batch_sz % numthreads == 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "batch size (-dm_landau_batch_size) %" PetscInt_FMT " mod #threads %" PetscInt_FMT " must equal zero", ctx->batch_sz, numthreads);
507: PetscCall(DMCompositeGetNumberDM(pack, &nDMs)); // number of vertices * number of grids
508: PetscCall(PetscInfo(pack, "Have %" PetscInt_FMT " total grids, with %" PetscInt_FMT " Landau local batched and %" PetscInt_FMT " global items (vertices) %d DMs\n", ctx->num_grids, ctx->batch_sz, num_vertices, (int)nDMs));
509: PetscCall(PetscMalloc(sizeof(*globXArray) * nDMs, &globXArray));
510: PetscCall(PetscMalloc(sizeof(*globMpArray) * nDMs, &globMpArray));
511: PetscCall(PetscMalloc(sizeof(*globSwarmArray) * nDMs, &globSwarmArray));
512: // print ctx
513: PetscCall(PetscNew(&printCtx));
514: PetscCall(TSSetApplicationContext(ts, printCtx));
515: printCtx->v_target = v_target;
516: printCtx->g_target = g_target;
517: printCtx->ctx = ctx;
518: printCtx->globSwarmArray = globSwarmArray;
519: printCtx->grid_dm = grid_dm;
520: printCtx->globMpArray = globMpArray;
521: printCtx->g_Mass = g_Mass;
522: printCtx->globXArray = globXArray;
523: printCtx->print_entropy = PETSC_FALSE;
524: PetscOptionsBegin(PETSC_COMM_SELF, "", "Print Options", "DMPLEX");
525: PetscCall(PetscOptionsBool("-print_entropy", "Print entropy and moments at each time step", "ex30.c", printCtx->print_entropy, &printCtx->print_entropy, NULL));
526: PetscOptionsEnd();
527: // view
528: PetscCall(DMViewFromOptions(ctx->plex[g_target], NULL, "-ex30_dm_view"));
529: if (ctx->num_grids > g_target + 1) { PetscCall(DMViewFromOptions(ctx->plex[g_target + 1], NULL, "-ex30_dm_view2")); }
530: // create mesh mass matrices
531: PetscCall(VecZeroEntries(X));
532: PetscCall(DMCompositeGetAccessArray(pack, X, nDMs, NULL, globXArray)); // just to duplicate
533: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
534: Vec subX = globXArray[LAND_PACK_IDX(0, grid)];
535: DM dm = ctx->plex[grid];
536: PetscSection s;
537: grid_dm[grid] = dm;
538: PetscCall(DMCreateMassMatrix(dm, dm, &g_Mass[grid]));
539: //
540: PetscCall(DMGetLocalSection(dm, &s));
541: PetscCall(DMPlexCreateClosureIndex(dm, s));
542: for (PetscInt tid = 0; tid < numthreads; tid++) {
543: PC pc;
544: PetscCall(VecDuplicate(subX, &t_fhat[grid][tid]));
545: PetscCall(KSPCreate(PETSC_COMM_SELF, &t_ksp[grid][tid]));
546: PetscCall(KSPSetType(t_ksp[grid][tid], KSPCG));
547: PetscCall(KSPGetPC(t_ksp[grid][tid], &pc));
548: PetscCall(PCSetType(pc, PCJACOBI));
549: PetscCall(KSPSetOptionsPrefix(t_ksp[grid][tid], "ptof_"));
550: PetscCall(KSPSetOperators(t_ksp[grid][tid], g_Mass[grid], g_Mass[grid]));
551: PetscCall(KSPSetFromOptions(t_ksp[grid][tid]));
552: }
553: }
554: PetscCall(DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, globXArray));
555: PetscCall(TSGetTimeStep(ts, &dt_init)); // we could have an adaptive time stepper
556: // loop over all vertices in chucks that are batched for TSSolve
557: for (PetscInt i = 0; i < 5; i++) moments_0[i] = moments_1a[i] = moments_1b[i] = 0;
558: for (PetscInt global_vertex_id_0 = 0; global_vertex_id_0 < num_vertices; global_vertex_id_0 += ctx->batch_sz, shift /= 2) { // outer vertex loop
559: PetscCall(TSSetTime(ts, 0));
560: PetscCall(TSSetStepNumber(ts, 0));
561: PetscCall(TSSetTimeStep(ts, dt_init));
562: PetscCall(DMCompositeGetAccessArray(pack, X, nDMs, NULL, globXArray));
563: printCtx->global_vertex_id_0 = global_vertex_id_0;
564: if (v_target >= global_vertex_id_0 && v_target < global_vertex_id_0 + ctx->batch_sz) {
565: PetscCall(PetscObjectSetName((PetscObject)globXArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target)], "rho"));
566: printCtx->print = PETSC_TRUE;
567: } else printCtx->print = PETSC_FALSE;
568: // create fake particles in batches with threads
569: for (PetscInt v_id_0 = 0; v_id_0 < ctx->batch_sz; v_id_0 += numthreads) {
570: PetscReal *xx_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS], *yy_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS], *zz_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS], *wp_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS] /* , radiuses[80000] */;
571: PetscInt Np_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS];
572: // make particles
573: for (PetscInt tid = 0; tid < numthreads; tid++) {
574: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
575: if (glb_v_id < num_vertices) { // the ragged edge (in last batch)
576: PetscInt Npp0 = a_Np + (glb_v_id % (a_Np / 10 + 1)), nTargetP[LANDAU_MAX_GRIDS]; // n of particels in each dim with load imbalance
577: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
578: // for (PetscInt sp = ctx->species_offset[grid], i0 = 0; sp < ctx->species_offset[grid + 1]; sp++, i0++) {
579: const PetscReal kT_m = ctx->k * ctx->thermal_temps[ctx->species_offset[grid]] / ctx->masses[ctx->species_offset[grid]] / (ctx->v_0 * ctx->v_0); /* theta = 2kT/mc^2 per species */
580: PetscReal lo[3] = {-ctx->radius[grid], -ctx->radius[grid], -ctx->radius[grid]}, hi[3] = {ctx->radius[grid], ctx->radius[grid], ctx->radius[grid]}, hp[3], vole; // would be nice to get box from DM
581: PetscInt Npi = Npp0, Npj = 2 * Npp0, Npk = 1;
582: PetscRandom rand;
583: PetscReal sigma = ctx->thermal_speed[grid] / ctx->thermal_speed[0], p2_shift = grid == 0 ? shift : -shift; // symmetric shift of e vs ions
584: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
585: PetscCall(PetscRandomSetInterval(rand, 0., 1.));
586: PetscCall(PetscRandomSetFromOptions(rand));
587: if (dim == 2) lo[0] = 0; // Landau coordinate (r,z)
588: else Npi = Npj = Npk = Npp0;
589: // User: use glb_v_id to index into your data
590: const PetscInt NNreal = Npi * Npj * Npk, NN = NNreal + (dim == 2 ? 3 : 6); // make room for bounding box
591: Np_t[grid][tid] = NN;
592: if (glb_v_id == v_target) nTargetP[grid] = NN;
593: PetscCall(PetscMalloc4(NN, &xx_t[grid][tid], NN, &yy_t[grid][tid], NN, &wp_t[grid][tid], dim == 2 ? 1 : NN, &zz_t[grid][tid]));
594: hp[0] = (hi[0] - lo[0]) / Npi;
595: hp[1] = (hi[1] - lo[1]) / Npj;
596: hp[2] = (hi[2] - lo[2]) / Npk;
597: if (dim == 2) hp[2] = 1;
598: PetscCall(PetscInfo(pack, " lo = %14.7e, hi = %14.7e; hp = %14.7e, %14.7e; kT_m = %g; \n", (double)lo[1], (double)hi[1], (double)hp[0], (double)hp[1], (double)kT_m)); // temp
599: vole = hp[0] * hp[1] * hp[2] * ctx->n[grid]; // fix for multi-species
600: PetscCall(PetscInfo(pack, "Vertex %" PetscInt_FMT ", grid %" PetscInt_FMT " with %" PetscInt_FMT " particles (diagnostic target = %" PetscInt_FMT ")\n", glb_v_id, grid, NN, v_target));
601: for (PetscInt pj = 0, pp = 0; pj < Npj; pj++) {
602: for (PetscInt pk = 0; pk < Npk; pk++) {
603: for (PetscInt pi = 0; pi < Npi; pi++, pp++) {
604: PetscReal p_shift = p2_shift;
605: wp_t[grid][tid][pp] = 0;
606: if (use_uniform_particle_grid) {
607: xx_t[grid][tid][pp] = lo[0] + hp[0] / 2.0 + pi * hp[0];
608: yy_t[grid][tid][pp] = lo[1] + hp[1] / 2.0 + pj * hp[1];
609: if (dim == 3) zz_t[grid][tid][pp] = lo[2] + hp[2] / 2.0 + pk * hp[2];
610: PetscReal x[] = {xx_t[grid][tid][pp], yy_t[grid][tid][pp], dim == 2 ? 0 : zz_t[grid][tid][pp]};
611: p_shift *= ctx->thermal_speed[grid] / ctx->v_0;
612: if (ctx->sphere && PetscSqrtReal(PetscSqr(xx_t[grid][tid][pp]) + PetscSqr(yy_t[grid][tid][pp])) > 0.92 * hi[0]) {
613: wp_t[grid][tid][pp] = 0;
614: } else {
615: maxwellian(dim, x, kT_m, vole, p_shift, &wp_t[grid][tid][pp]);
616: if (ctx->num_grids == 1 && shift != 0) { // bi-maxwellian, electron plasma
617: maxwellian(dim, x, kT_m, vole, -p_shift, &wp_t[grid][tid][pp]); // symmetric shift of electron plasma
618: }
619: }
620: } else {
621: PetscReal u1, u2;
622: do {
623: do {
624: PetscCall(PetscRandomGetValueReal(rand, &u1));
625: } while (u1 == 0);
626: PetscCall(PetscRandomGetValueReal(rand, &u2));
627: //compute z0 and z1
628: PetscReal mag = sigma * PetscSqrtReal(-2.0 * PetscLogReal(u1)); // is this the same scale grid Maxwellian? t_therm = sigma
629: xx_t[grid][tid][pp] = mag * PetscCosReal(2.0 * PETSC_PI * u2);
630: yy_t[grid][tid][pp] = mag * PetscSinReal(2.0 * PETSC_PI * u2);
631: if (dim == 2 && xx_t[grid][tid][pp] < lo[0]) xx_t[grid][tid][pp] = -xx_t[grid][tid][pp];
632: if (dim == 3) zz_t[grid][tid][pp] = lo[2] + hp[2] / 2.0 + pk * hp[2];
633: if (!ctx->sphere) {
634: if (dim == 2 && xx_t[grid][tid][pp] < 0) xx_t[grid][tid][pp] = -xx_t[grid][tid][pp]; // ???
635: else if (dim == 3) {
636: while (zz_t[grid][tid][pp] >= hi[2] || zz_t[grid][tid][pp] <= lo[2]) zz_t[grid][tid][pp] *= .9;
637: }
638: while (xx_t[grid][tid][pp] >= hi[0] || xx_t[grid][tid][pp] <= lo[0]) xx_t[grid][tid][pp] *= .9;
639: while (yy_t[grid][tid][pp] >= hi[1] || yy_t[grid][tid][pp] <= lo[1]) yy_t[grid][tid][pp] *= .9;
640: } else { // 2D
641: //if (glb_v_id == v_target && pp < 80000) radiuses[pp] = PetscSqrtReal(PetscSqr(xx_t[grid][tid][pp]) + PetscSqr(yy_t[grid][tid][pp]));
642: while (PetscSqrtReal(PetscSqr(xx_t[grid][tid][pp]) + PetscSqr(yy_t[grid][tid][pp])) > 0.92 * hi[0]) { // safety factor for facets of sphere
643: xx_t[grid][tid][pp] *= .9;
644: yy_t[grid][tid][pp] *= .9;
645: }
646: }
647: if (ctx->num_grids == 1 && pp % 2 == 0) p_shift = 0; // one species, split bi-max
648: p_shift *= ctx->thermal_speed[grid] / ctx->v_0;
649: if (dim == 3) zz_t[grid][tid][pp] += p_shift;
650: else yy_t[grid][tid][pp] += p_shift;
651: wp_t[grid][tid][pp] += ctx->n[grid] / NNreal * PetscSqrtReal(ctx->masses[ctx->species_offset[grid]] / ctx->masses[0]);
652: if (p_shift <= 0) break; // add bi-max for electron plasma only
653: p_shift = -p_shift;
654: } while (ctx->num_grids == 1); // add bi-max for electron plasma only
655: }
656: {
657: if (glb_v_id == v_target) {
658: PetscReal x[] = {xx_t[grid][tid][pp], yy_t[grid][tid][pp], dim == 2 ? 0 : zz_t[grid][tid][pp]};
659: PetscReal v2 = 0, fact = dim == 2 ? 2.0 * PETSC_PI * x[0] : 1, w = fact * wp_t[grid][tid][pp] * ctx->n_0 * ctx->masses[ctx->species_offset[grid]];
660: for (PetscInt i = 0; i < dim; ++i) v2 += PetscSqr(x[i]);
661: moments_0[0] += w; // not thread safe
662: moments_0[1] += w * ctx->v_0 * x[1]; // z-momentum
663: moments_0[2] += w * 0.5 * ctx->v_0 * ctx->v_0 * v2;
664: }
665: }
666: }
667: }
668: }
669: if (dim == 2) { // fix bounding box
670: PetscInt pp = NNreal;
671: wp_t[grid][tid][pp] = 0;
672: xx_t[grid][tid][pp] = 1.e-7;
673: yy_t[grid][tid][pp++] = hi[1] - 5.e-7;
674: wp_t[grid][tid][pp] = 0;
675: xx_t[grid][tid][pp] = hi[0] - 5.e-7;
676: yy_t[grid][tid][pp++] = 0;
677: wp_t[grid][tid][pp] = 0;
678: xx_t[grid][tid][pp] = 1.e-7;
679: yy_t[grid][tid][pp++] = lo[1] + 5.e-7;
680: } else {
681: const PetscInt p0 = NNreal;
682: for (PetscInt pj = 0; pj < 6; pj++) { xx_t[grid][tid][p0 + pj] = yy_t[grid][tid][p0 + pj] = zz_t[grid][tid][p0 + pj] = wp_t[grid][tid][p0 + pj] = 0; }
683: xx_t[grid][tid][p0 + 0] = lo[0];
684: xx_t[grid][tid][p0 + 1] = hi[0];
685: yy_t[grid][tid][p0 + 2] = lo[1];
686: yy_t[grid][tid][p0 + 3] = hi[1];
687: zz_t[grid][tid][p0 + 4] = lo[2];
688: zz_t[grid][tid][p0 + 5] = hi[2];
689: }
690: PetscCall(PetscRandomDestroy(&rand));
691: }
692: // entropy init, need global n
693: if (glb_v_id == v_target) {
694: const PetscReal N_inv = 1 / moments_0[0];
695: PetscCall(PetscInfo(pack, "Target %" PetscInt_FMT " with %" PetscInt_FMT " particels\n", glb_v_id, nTargetP[0]));
696: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
697: const PetscInt NN = nTargetP[grid];
698: for (PetscInt pp = 0; pp < NN; pp++) {
699: const PetscReal fact = dim == 2 ? 2.0 * PETSC_PI * xx_t[grid][tid][pp] : 1, w = fact * ctx->n_0 * ctx->masses[ctx->species_offset[grid]] * wp_t[grid][tid][pp], ww = w * N_inv;
700: if (w > PETSC_REAL_MIN) {
701: moments_0[3] -= ww * PetscLogReal(ww);
702: PetscCheck(ww < 1 - PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "ww (%g) > 1", (double)ww);
703: } else moments_0[4] -= w;
704: }
705: } // grid
706: } // target
707: } // active
708: } // threads
709: /* Create particle swarm */
710: for (PetscInt tid = 0; tid < numthreads; tid++) {
711: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
712: if (glb_v_id < num_vertices) { // the ragged edge of the last batch
713: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
714: PetscSection section;
715: PetscInt Nf;
716: DM dm = grid_dm[grid];
717: PetscCall(DMGetLocalSection(dm, §ion));
718: PetscCall(PetscSectionGetNumFields(section, &Nf));
719: PetscCheck(Nf == 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Only one species per grid supported -- todo");
720: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
721: PetscCall(PetscInfo(pack, "call createSwarm [%" PetscInt_FMT ".%" PetscInt_FMT "] local block index %" PetscInt_FMT "\n", v_id, grid, LAND_PACK_IDX(v_id, grid)));
722: PetscCall(createSwarm(dm, dim, &globSwarmArray[LAND_PACK_IDX(v_id, grid)]));
723: }
724: } // active
725: } // threads
726: PetscCheck(ierr != 9999, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Only support one species per grid");
727: // make globMpArray
728: PetscPragmaOMP(parallel for)
729: for (PetscInt tid = 0; tid < numthreads; tid++) {
730: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
731: if (glb_v_id < num_vertices) {
732: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
733: // for (PetscInt sp = ctx->species_offset[grid], i0 = 0; sp < ctx->species_offset[grid + 1]; sp++, i0++) -- loop over species for Nf > 1 -- TODO
734: PetscErrorCode ierr_t;
735: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
736: ierr_t = PetscInfo(pack, "makeSwarm %" PetscInt_FMT ".%" PetscInt_FMT ") for block %" PetscInt_FMT "\n", v_id, grid, LAND_PACK_IDX(v_id, grid));
737: ierr_t = makeSwarm(sw, dim, Np_t[grid][tid], xx_t[grid][tid], yy_t[grid][tid], zz_t[grid][tid]);
738: if (ierr_t) ierr = ierr_t;
739: }
740: }
741: }
742: for (PetscInt tid = 0; tid < numthreads; tid++) {
743: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
744: if (glb_v_id < num_vertices) {
745: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
746: DM dm = grid_dm[grid];
747: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
748: PetscCall(PetscInfo(pack, "createMp %" PetscInt_FMT ".%" PetscInt_FMT ") for block %" PetscInt_FMT "\n", v_id, grid, LAND_PACK_IDX(v_id, grid)));
749: PetscCall(createMp(dm, sw, &globMpArray[LAND_PACK_IDX(v_id, grid)]));
750: PetscCall(MatViewFromOptions(globMpArray[LAND_PACK_IDX(v_id, grid)], NULL, "-mp_mat_view"));
751: }
752: }
753: }
754: // p --> g: set X
755: // PetscPragmaOMP(parallel for)
756: for (PetscInt tid = 0; tid < numthreads; tid++) {
757: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
758: if (glb_v_id < num_vertices) {
759: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
760: PetscErrorCode ierr_t;
761: DM dm = grid_dm[grid];
762: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
763: Vec subX = globXArray[LAND_PACK_IDX(v_id, grid)], work = t_fhat[grid][tid];
764: ierr_t = PetscInfo(pack, "particlesToGrid %" PetscInt_FMT ".%" PetscInt_FMT ") for block %" PetscInt_FMT "\n", v_id, grid, LAND_PACK_IDX(v_id, grid));
765: ierr_t = particlesToGrid(dm, sw, tid, dim, wp_t[grid][tid], subX, globMpArray[LAND_PACK_IDX(v_id, grid)]);
766: if (ierr_t) ierr = ierr_t;
767: // u = M^_1 f_w
768: ierr_t = VecCopy(subX, work);
769: ierr_t = KSPSolve(t_ksp[grid][tid], work, subX);
770: if (ierr_t) ierr = ierr_t;
771: }
772: }
773: }
774: PetscCheck(!ierr, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in OMP loop. ierr = %d", (int)ierr);
775: /* Cleanup */
776: for (PetscInt tid = 0; tid < numthreads; tid++) {
777: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
778: if (glb_v_id < num_vertices) {
779: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
780: PetscCall(PetscFree4(xx_t[grid][tid], yy_t[grid][tid], wp_t[grid][tid], zz_t[grid][tid]));
781: }
782: } // active
783: } // threads
784: } // (fake) particle loop
785: // standard view of initial conditions
786: if (v_target >= global_vertex_id_0 && v_target < global_vertex_id_0 + ctx->batch_sz) {
787: PetscCall(DMSetOutputSequenceNumber(ctx->plex[g_target], 0, 0.0));
788: PetscCall(VecViewFromOptions(globXArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target)], NULL, "-ex30_vec_view"));
789: if (ctx->num_grids > g_target + 1) {
790: PetscCall(DMSetOutputSequenceNumber(ctx->plex[g_target + 1], 0, 0.0));
791: PetscCall(VecViewFromOptions(globXArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target + 1)], NULL, "-ex30_vec_view2"));
792: }
793: PetscCall(MatViewFromOptions(globMpArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target)], NULL, "-ex30_mass_mat_view"));
794: PetscCall(DMViewFromOptions(globSwarmArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target)], NULL, "-ex30_sw_view"));
795: PetscCall(DMSwarmViewXDMF(globSwarmArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target)], "initial_swarm.xmf")); // writes a file by default!!!
796: }
797: // coarse graining moments_1a, bring f back from grid before advance
798: if (v_target >= global_vertex_id_0 && v_target < global_vertex_id_0 + ctx->batch_sz && printCtx->print_entropy) {
799: const PetscInt v_id = v_target % ctx->batch_sz;
800: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
801: PetscDataType dtype;
802: PetscReal *wp, *coords;
803: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
804: Vec work, subX = globXArray[LAND_PACK_IDX(v_id, grid)];
805: PetscInt bs, NN;
806: // C-G moments
807: PetscCall(VecDuplicate(subX, &work));
808: PetscCall(gridToParticles(grid_dm[grid], sw, subX, work, globMpArray[LAND_PACK_IDX(v_id, grid)], g_Mass[grid]));
809: PetscCall(VecDestroy(&work));
810: // moments
811: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
812: PetscCall(DMSwarmGetLocalSize(sw, &NN));
813: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wp));
814: for (PetscInt pp = 0; pp < NN; pp++) {
815: PetscReal v2 = 0, fact = (dim == 2) ? 2.0 * PETSC_PI * coords[pp * dim + 0] : 1, w = fact * wp[pp] * ctx->n_0 * ctx->masses[ctx->species_offset[grid]];
816: for (PetscInt i = 0; i < dim; ++i) v2 += PetscSqr(coords[pp * dim + i]);
817: moments_1a[0] += w;
818: moments_1a[1] += w * ctx->v_0 * coords[pp * dim + 1]; // z-momentum
819: moments_1a[2] += w * 0.5 * ctx->v_0 * ctx->v_0 * v2;
820: }
821: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
822: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wp));
823: }
824: // entropy
825: const PetscReal N_inv = 1 / moments_1a[0];
826: PetscCall(PetscInfo(pack, "Entropy batch %" PetscInt_FMT " of %" PetscInt_FMT ", n = %g\n", v_target, num_vertices, (double)(1 / N_inv)));
827: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
828: PetscDataType dtype;
829: PetscReal *wp, *coords;
830: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
831: PetscInt bs, NN;
832: PetscCall(DMSwarmGetLocalSize(sw, &NN));
833: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wp));
834: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
835: for (PetscInt pp = 0; pp < NN; pp++) {
836: PetscReal fact = (dim == 2) ? 2.0 * PETSC_PI * coords[pp * dim + 0] : 1, w = fact * wp[pp] * ctx->n_0 * ctx->masses[ctx->species_offset[grid]], ww = w * N_inv;
837: if (w > PETSC_REAL_MIN) {
838: moments_1a[3] -= ww * PetscLogReal(ww);
839: PetscCheck(ww < 1 - PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "ww (%g) > 1", (double)ww);
840: } else moments_1a[4] -= w;
841: }
842: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wp));
843: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
844: }
845: }
846: // restore vector
847: PetscCall(DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, globXArray));
848: // view initial grid
849: if (v_target >= global_vertex_id_0 && v_target < global_vertex_id_0 + ctx->batch_sz) { PetscCall(DMPlexLandauPrintNorms(X, 0)); }
850: // advance
851: PetscCall(TSSetSolution(ts, X));
852: PetscCall(PetscInfo(pack, "Advance vertex %" PetscInt_FMT " to %" PetscInt_FMT "\n", global_vertex_id_0, global_vertex_id_0 + ctx->batch_sz));
853: PetscCall(TSSetPostStep(ts, PostStep));
854: PetscCall(PostStep(ts));
855: PetscCall(TSSolve(ts, X));
856: // view
857: PetscCall(DMCompositeGetAccessArray(pack, X, nDMs, NULL, globXArray));
858: if (v_target >= global_vertex_id_0 && v_target < global_vertex_id_0 + ctx->batch_sz) {
859: /* Visualize original particle field */
860: DM sw = globSwarmArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target)];
861: Vec f;
862: PetscCall(DMSetOutputSequenceNumber(sw, 0, 0.0));
863: PetscCall(DMViewFromOptions(grid_dm[g_target], NULL, "-weights_dm_view"));
864: PetscCall(DMViewFromOptions(sw, NULL, "-weights_sw_view"));
865: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
866: PetscCall(PetscObjectSetName((PetscObject)f, "weights"));
867: PetscCall(VecViewFromOptions(f, NULL, "-weights_vec_view"));
868: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
869: //
870: PetscCall(DMPlexLandauPrintNorms(X, 1));
871: }
872: if (!use_uniform_particle_grid) { // resample to uniform grid
873: for (PetscInt v_id_0 = 0; v_id_0 < ctx->batch_sz; v_id_0 += numthreads) {
874: PetscReal *xx_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS], *yy_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS], *zz_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS], *wp_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS];
875: PetscInt Np_t[LANDAU_MAX_GRIDS][EX30_MAX_NUM_THRDS];
876: for (PetscInt tid = 0; tid < numthreads; tid++) {
877: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
878: if (glb_v_id < num_vertices) {
879: // create uniform grid w/o weights & smaller
880: PetscInt Npp0 = (a_Np + (glb_v_id % (a_Np / 10 + 1))) / 2, Nv; // 1/2 of uniform particle grid size
881: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
882: // for (PetscInt sp = ctx->species_offset[grid], i0 = 0; sp < ctx->species_offset[grid + 1]; sp++, i0++)
883: PetscReal lo[3] = {-ctx->radius[grid], -ctx->radius[grid], -ctx->radius[grid]}, hi[3] = {ctx->radius[grid], ctx->radius[grid], ctx->radius[grid]}, hp[3];
884: PetscInt Npi = Npp0, Npj = 2 * Npp0, Npk = 1, NN;
885: // delete old particles and particle mass matrix
886: PetscCall(DMDestroy(&globSwarmArray[LAND_PACK_IDX(v_id, grid)]));
887: PetscCall(MatDestroy(&globMpArray[LAND_PACK_IDX(v_id, grid)]));
888: // create fake particles in batches with threads
889: PetscCall(MatGetLocalSize(g_Mass[grid], &Nv, NULL));
890: if (dim == 2) lo[0] = 0;
891: else Npi = Npj = Npk = Npp0;
892: NN = Npi * Npj * Npk + (dim == 2 ? 3 : 6); // make a regular grid of particles Npp x Npp
893: while (Npi * Npj * Npk < Nv) { // make stable - no LS
894: Npi++;
895: Npj++;
896: Npk++;
897: NN = Npi * Npj * Npk + (dim == 2 ? 3 : 6);
898: }
899: Np_t[grid][tid] = NN;
900: PetscCall(PetscMalloc4(NN, &xx_t[grid][tid], NN, &yy_t[grid][tid], NN, &wp_t[grid][tid], dim == 2 ? 1 : NN, &zz_t[grid][tid]));
901: hp[0] = (hi[0] - lo[0]) / Npi;
902: hp[1] = (hi[1] - lo[1]) / Npj;
903: hp[2] = (hi[2] - lo[2]) / Npk;
904: if (dim == 2) hp[2] = 1;
905: PetscCall(PetscInfo(pack, "Resampling %d particles, %d vertices\n", (int)NN, (int)Nv)); // temp
906: for (PetscInt pj = 0, pp = 0; pj < Npj; pj++) {
907: for (PetscInt pk = 0; pk < Npk; pk++) {
908: for (PetscInt pi = 0; pi < Npi; pi++, pp++) {
909: wp_t[grid][tid][pp] = 0;
910: xx_t[grid][tid][pp] = lo[0] + hp[0] / 2.0 + pi * hp[0];
911: yy_t[grid][tid][pp] = lo[1] + hp[1] / 2.0 + pj * hp[1];
912: if (dim == 3) zz_t[grid][tid][pp] = lo[2] + hp[2] / 2.0 + pk * hp[2];
913: }
914: }
915: }
916: if (dim == 2) { // fix bounding box
917: PetscInt pp = NN - 3;
918: wp_t[grid][tid][pp] = 0;
919: xx_t[grid][tid][pp] = 1.e-7;
920: yy_t[grid][tid][pp++] = hi[1] - 5.e-7;
921: wp_t[grid][tid][pp] = 0;
922: xx_t[grid][tid][pp] = hi[0] - 5.e-7;
923: yy_t[grid][tid][pp++] = 0;
924: wp_t[grid][tid][pp] = 0;
925: xx_t[grid][tid][pp] = 1.e-7;
926: yy_t[grid][tid][pp++] = lo[1] + 5.e-7;
927: } else {
928: const PetscInt p0 = NN - 6;
929: for (PetscInt pj = 0; pj < 6; pj++) { xx_t[grid][tid][p0 + pj] = yy_t[grid][tid][p0 + pj] = zz_t[grid][tid][p0 + pj] = wp_t[grid][tid][p0 + pj] = 0; }
930: xx_t[grid][tid][p0 + 0] = lo[0];
931: xx_t[grid][tid][p0 + 1] = hi[0];
932: yy_t[grid][tid][p0 + 2] = lo[1];
933: yy_t[grid][tid][p0 + 3] = hi[1];
934: zz_t[grid][tid][p0 + 4] = lo[2];
935: zz_t[grid][tid][p0 + 5] = hi[2];
936: }
937: }
938: } // active
939: } // threads
940: /* Create particle swarm */
941: for (PetscInt tid = 0; tid < numthreads; tid++) {
942: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
943: if (glb_v_id < num_vertices) { // the ragged edge of the last batch
944: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
945: // for (PetscInt sp = ctx->species_offset[grid], i0 = 0; sp < ctx->species_offset[grid + 1]; sp++, i0++) -- loop over species for Nf > 1 -- TODO
946: PetscErrorCode ierr_t;
947: PetscSection section;
948: PetscInt Nf;
949: DM dm = grid_dm[grid];
950: ierr_t = DMGetLocalSection(dm, §ion);
951: ierr_t = PetscSectionGetNumFields(section, &Nf);
952: if (Nf != 1) ierr_t = (PetscErrorCode)9999;
953: else {
954: ierr_t = DMViewFromOptions(dm, NULL, "-dm_view");
955: ierr_t = PetscInfo(pack, "call createSwarm [%" PetscInt_FMT ".%" PetscInt_FMT "] local block index %" PetscInt_FMT "\n", v_id, grid, LAND_PACK_IDX(v_id, grid));
956: ierr_t = createSwarm(dm, dim, &globSwarmArray[LAND_PACK_IDX(v_id, grid)]);
957: }
958: if (ierr_t) ierr = ierr_t;
959: }
960: } // active
961: } // threads
962: PetscCheck(ierr != 9999, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Only support one species per grid");
963: PetscCheck(!ierr, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Error in OMP loop. ierr = %d", (int)ierr);
964: // make globMpArray
965: PetscPragmaOMP(parallel for)
966: for (PetscInt tid = 0; tid < numthreads; tid++) {
967: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
968: if (glb_v_id < num_vertices) {
969: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
970: // for (PetscInt sp = ctx->species_offset[grid], i0 = 0; sp < ctx->species_offset[grid + 1]; sp++, i0++) -- loop over species for Nf > 1 -- TODO
971: PetscErrorCode ierr_t;
972: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
973: ierr_t = PetscInfo(pack, "makeSwarm %" PetscInt_FMT ".%" PetscInt_FMT ") for block %" PetscInt_FMT "\n", v_id, grid, LAND_PACK_IDX(v_id, grid));
974: ierr_t = makeSwarm(sw, dim, Np_t[grid][tid], xx_t[grid][tid], yy_t[grid][tid], zz_t[grid][tid]);
975: if (ierr_t) ierr = ierr_t;
976: }
977: } // active
978: } // threads
979: // create particle mass matrices
980: //PetscPragmaOMP(parallel for)
981: for (PetscInt tid = 0; tid < numthreads; tid++) {
982: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
983: if (glb_v_id < num_vertices) {
984: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
985: PetscErrorCode ierr_t;
986: DM dm = grid_dm[grid];
987: DM sw = globSwarmArray[LAND_PACK_IDX(v_id, grid)];
988: ierr_t = PetscInfo(pack, "createMp %" PetscInt_FMT ".%" PetscInt_FMT ") for block %" PetscInt_FMT "\n", v_id, grid, LAND_PACK_IDX(v_id, grid));
989: ierr_t = createMp(dm, sw, &globMpArray[LAND_PACK_IDX(v_id, grid)]);
990: if (ierr_t) ierr = ierr_t;
991: }
992: } // active
993: } // threads
994: PetscCheck(!ierr, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in OMP loop. ierr = %d", (int)ierr);
995: /* Cleanup */
996: for (PetscInt tid = 0; tid < numthreads; tid++) {
997: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
998: if (glb_v_id < num_vertices) {
999: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
1000: PetscCall(PetscFree4(xx_t[grid][tid], yy_t[grid][tid], wp_t[grid][tid], zz_t[grid][tid]));
1001: }
1002: } // active
1003: } // threads
1004: } // batch
1005: // view
1006: if (v_target >= global_vertex_id_0 && v_target < global_vertex_id_0 + ctx->batch_sz) {
1007: /* Visualize particle field */
1008: DM sw = globSwarmArray[LAND_PACK_IDX(v_target % ctx->batch_sz, g_target)];
1009: Vec f;
1010: PetscCall(DMSetOutputSequenceNumber(sw, 0, 0.0));
1011: PetscCall(DMViewFromOptions(sw, NULL, "-resampled_weights_sw_view"));
1012: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
1013: PetscCall(PetscObjectSetName((PetscObject)f, "resampled_weights"));
1014: PetscCall(VecViewFromOptions(f, NULL, "-resampled_weights_vec_view"));
1015: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
1016: PetscCall(DMSwarmViewXDMF(sw, "resampled.xmf"));
1017: }
1018: } // !uniform
1019: // particles to grid, compute moments and entropy, for target vertex only
1020: if (v_target >= global_vertex_id_0 && v_target < global_vertex_id_0 + ctx->batch_sz && printCtx->print_entropy) {
1021: PetscReal energy_error_rel;
1022: PetscCall(gridToParticles_private(grid_dm, globSwarmArray, dim, v_target, numthreads, num_vertices, global_vertex_id_0, globMpArray, g_Mass, t_fhat, moments_1b, globXArray, ctx));
1023: energy_error_rel = PetscAbsReal(moments_1b[2] - moments_0[2]) / moments_0[2];
1024: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Particle Moments:\t number density momentum (par) energy entropy negative weights : # OMP threads %g\n", (double)numthreads));
1025: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\tInitial: %18.12e %19.12e %18.12e %18.12e %g %%\n", (double)moments_0[0], (double)moments_0[1], (double)moments_0[2], (double)moments_0[3], 100 * (double)(moments_0[4] / moments_0[0])));
1026: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\tCoarse-graining: %18.12e %19.12e %18.12e %18.12e %g %%\n", (double)moments_1a[0], (double)moments_1a[1], (double)moments_1a[2], (double)moments_1a[3], 100 * (double)(moments_1a[4] / moments_0[0])));
1027: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\tLandau: %18.12e %19.12e %18.12e %18.12e %g %%\n", (double)moments_1b[0], (double)moments_1b[1], (double)moments_1b[2], (double)moments_1b[3], 100 * (double)(moments_1b[4] / moments_0[0])));
1028: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coarse-graining entropy generation = %e ; Landau entropy generation = %e\n", (double)(moments_1a[3] - moments_0[3]), (double)(moments_1b[3] - moments_0[3])));
1029: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(relative) energy conservation: Coarse-graining = %e, Landau = %e (%g %d)\n", (double)(moments_1a[2] - moments_0[2]) / (double)moments_0[2], (double)energy_error_rel, (double)PetscLog10Real(energy_error_rel), (int)(PetscLog10Real(energy_error_rel) + .5)));
1030: }
1031: // restore vector
1032: PetscCall(DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, globXArray));
1033: // cleanup
1034: for (PetscInt v_id_0 = 0; v_id_0 < ctx->batch_sz; v_id_0 += numthreads) {
1035: for (PetscInt tid = 0; tid < numthreads; tid++) {
1036: const PetscInt v_id = v_id_0 + tid, glb_v_id = global_vertex_id_0 + v_id;
1037: if (glb_v_id < num_vertices) {
1038: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
1039: PetscCall(DMDestroy(&globSwarmArray[LAND_PACK_IDX(v_id, grid)]));
1040: PetscCall(MatDestroy(&globMpArray[LAND_PACK_IDX(v_id, grid)]));
1041: }
1042: }
1043: }
1044: }
1045: } // user batch, not used
1046: /* Cleanup */
1047: PetscCall(PetscFree(globXArray));
1048: PetscCall(PetscFree(globSwarmArray));
1049: PetscCall(PetscFree(globMpArray));
1050: PetscCall(PetscFree(printCtx));
1051: // clean up mass matrices
1052: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // add same particels for all grids
1053: PetscCall(MatDestroy(&g_Mass[grid]));
1054: for (PetscInt tid = 0; tid < numthreads; tid++) {
1055: PetscCall(VecDestroy(&t_fhat[grid][tid]));
1056: PetscCall(KSPDestroy(&t_ksp[grid][tid]));
1057: }
1058: }
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: int main(int argc, char **argv)
1063: {
1064: DM pack;
1065: Vec X;
1066: PetscInt dim = 2, num_vertices = 1, Np = 10, v_target = 0, g_target = 0;
1067: TS ts;
1068: Mat J;
1069: LandauCtx *ctx;
1070: PetscReal shift = 0;
1071: PetscBool use_uniform_particle_grid = PETSC_TRUE;
1073: PetscFunctionBeginUser;
1074: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1075: // process args
1076: PetscOptionsBegin(PETSC_COMM_SELF, "", "Collision Options", "DMPLEX");
1077: PetscCall(PetscOptionsInt("-dim", "Velocity space dimension", "ex30.c", dim, &dim, NULL));
1078: PetscCall(PetscOptionsInt("-number_spatial_vertices", "Number of user spatial vertices to be batched for Landau", "ex30.c", num_vertices, &num_vertices, NULL));
1079: PetscCall(PetscOptionsInt("-number_particles_per_dimension", "Number of particles per grid, with slight modification per spatial vertex, in each dimension of base Cartesian grid", "ex30.c", Np, &Np, NULL));
1080: PetscCall(PetscOptionsBool("-use_uniform_particle_grid", "Use uniform particle grid", "ex30.c", use_uniform_particle_grid, &use_uniform_particle_grid, NULL));
1081: PetscCall(PetscOptionsInt("-vertex_view_target", "Global vertex for diagnostics", "ex30.c", v_target, &v_target, NULL));
1082: PetscCall(PetscOptionsReal("-e_shift", "Bi-Maxwellian shift", "ex30.c", shift, &shift, NULL));
1083: PetscCheck(v_target < num_vertices, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Batch to view %" PetscInt_FMT " should be < number of vertices %" PetscInt_FMT, v_target, num_vertices);
1084: PetscCall(PetscOptionsInt("-grid_view_target", "Grid to view with diagnostics", "ex30.c", g_target, &g_target, NULL));
1085: PetscOptionsEnd();
1086: /* Create a mesh */
1087: PetscCall(DMPlexLandauCreateVelocitySpace(PETSC_COMM_SELF, dim, "", &X, &J, &pack));
1088: PetscCall(DMGetApplicationContext(pack, &ctx));
1089: PetscCall(DMSetUp(pack));
1090: PetscCall(DMSetOutputSequenceNumber(pack, 0, 0.0));
1091: PetscCheck(g_target < ctx->num_grids, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Grid to view %" PetscInt_FMT " should be < number of grids %" PetscInt_FMT, g_target, ctx->num_grids);
1092: PetscCheck(ctx->batch_view_idx == v_target % ctx->batch_sz, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Global view index %" PetscInt_FMT " mode batch size %" PetscInt_FMT " != ctx->batch_view_idx %" PetscInt_FMT, v_target, ctx->batch_sz, ctx->batch_view_idx);
1093: // PetscCheck(!use_uniform_particle_grid || !ctx->sphere, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Can not use -use_uniform_particle_grid and -dm_landau_sphere");
1094: /* Create timestepping solver context */
1095: PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
1096: PetscCall(TSSetDM(ts, pack));
1097: PetscCall(TSSetIFunction(ts, NULL, DMPlexLandauIFunction, NULL));
1098: PetscCall(TSSetIJacobian(ts, J, J, DMPlexLandauIJacobian, NULL));
1099: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1100: PetscCall(TSSetFromOptions(ts));
1101: PetscCall(PetscObjectSetName((PetscObject)X, "X"));
1102: // do particle advance
1103: PetscCall(go(ts, X, num_vertices, Np, dim, v_target, g_target, shift, use_uniform_particle_grid));
1104: PetscCall(MatZeroEntries(J)); // need to zero out so as to not reuse it in Landau's logic
1105: /* clean up */
1106: PetscCall(DMPlexLandauDestroyVelocitySpace(&pack));
1107: PetscCall(TSDestroy(&ts));
1108: PetscCall(VecDestroy(&X));
1109: PetscCall(PetscFinalize());
1110: return 0;
1111: }
1113: /*TEST
1115: build:
1116: requires: !complex
1118: testset:
1119: requires: double defined(PETSC_USE_DMLANDAU_2D)
1120: output_file: output/ex30_0.out
1121: args: -dim 2 -petscspace_degree 3 -dm_landau_num_species_grid 1,1,1 -dm_refine 1 -number_particles_per_dimension 20 \
1122: -dm_landau_batch_size 4 -number_spatial_vertices 6 -vertex_view_target 5 -grid_view_target 1 -dm_landau_batch_view_idx 1 \
1123: -dm_landau_n 1.000018,1,1e-6 -dm_landau_thermal_temps 2,1,1 -dm_landau_ion_masses 2,180 -dm_landau_ion_charges 1,18 \
1124: -ftop_ksp_rtol 1e-10 -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_factor_shift_type nonzero -ftop_sub_pc_type lu -ftop_ksp_error_if_not_converged \
1125: -ksp_type gmres -ksp_error_if_not_converged -dm_landau_verbose 4 -print_entropy \
1126: -ptof_ksp_type cg -ptof_pc_type jacobi -ptof_ksp_rtol 1e-12 -ptof_ksp_error_if_not_converged\
1127: -snes_converged_reason -snes_monitor -snes_rtol 1e-12 -snes_stol 1e-12 \
1128: -ts_dt 0.01 -ts_rtol 1e-1 -ts_exact_final_time stepover -ts_max_snes_failures -1 -ts_max_steps 1 -ts_monitor -ts_type beuler
1129: test:
1130: suffix: cpu
1131: args: -dm_landau_device_type cpu -pc_type jacobi
1132: test:
1133: suffix: kokkos
1134: # failed on Sunspot@ALCF with sycl
1135: requires: kokkos_kernels !openmp !sycl
1136: args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos -pc_type bjkokkos -pc_bjkokkos_ksp_type tfqmr -pc_bjkokkos_pc_type jacobi
1138: testset:
1139: requires: double !defined(PETSC_USE_DMLANDAU_2D)
1140: output_file: output/ex30_3d.out
1141: args: -dim 3 -petscspace_degree 2 -dm_landau_num_species_grid 1,1 -dm_refine 0 -number_particles_per_dimension 10 -dm_plex_hash_location \
1142: -dm_landau_batch_size 1 -number_spatial_vertices 1 -vertex_view_target 0 -grid_view_target 0 -dm_landau_batch_view_idx 0 \
1143: -dm_landau_n 1.000018,1 -dm_landau_thermal_temps 2,1 -dm_landau_ion_masses 2 -dm_landau_ion_charges 1 \
1144: -ftop_ksp_type cg -ftop_pc_type jacobi -ftop_ksp_rtol 1e-12 -ftop_ksp_error_if_not_converged -ksp_type preonly -pc_type lu -ksp_error_if_not_converged \
1145: -ptof_ksp_type cg -ptof_pc_type jacobi -ptof_ksp_rtol 1e-12 -ptof_ksp_error_if_not_converged \
1146: -snes_converged_reason -snes_monitor -snes_rtol 1e-12 -snes_stol 1e-12 \
1147: -ts_dt 0.1 -ts_exact_final_time stepover -ts_max_snes_failures -1 -ts_max_steps 1 -ts_monitor -ts_type beuler -print_entropy
1148: test:
1149: suffix: cpu_3d
1150: args: -dm_landau_device_type cpu
1151: test:
1152: suffix: kokkos_3d
1153: requires: kokkos_kernels !openmp
1154: args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos -pc_type bjkokkos -pc_bjkokkos_ksp_type tfqmr -pc_bjkokkos_pc_type jacobi
1156: test:
1157: suffix: conserve
1158: requires: !complex double defined(PETSC_USE_DMLANDAU_2D) !cuda
1159: args: -dm_landau_batch_size 4 -dm_refine 0 -dm_landau_num_species_grid 1 -dm_landau_thermal_temps 1 -petscspace_degree 3 -snes_converged_reason -ts_type beuler -ts_dt .1 \
1160: -ts_max_steps 1 -ksp_type preonly -ksp_error_if_not_converged -snes_rtol 1e-14 -snes_stol 1e-14 -dm_landau_device_type cpu -number_particles_per_dimension 20 \
1161: -ptof_ksp_type cg -ptof_pc_type jacobi -ptof_ksp_rtol 1e-14 -ptof_ksp_error_if_not_converged -pc_type lu -dm_landau_simplex 1 -use_uniform_particle_grid false -dm_landau_sphere -print_entropy -number_particles_per_dimension 50 -ftop_ksp_type cg -ftop_pc_type jacobi -ftop_ksp_rtol 1e-14
1163: TEST*/