Actual source code: ex2.c
1: static char help[] = "Runaway electron model with Landau collision operator\n\n";
3: #include <petscdmplex.h>
4: #include <petsclandau.h>
5: #include <petscts.h>
6: #include <petscds.h>
7: #include <petscdmcomposite.h>
8: #include <petsc/private/petscimpl.h>
10: #if defined(PETSC_HAVE_CUDA_NVTX)
11: #include <nvToolsExt.h>
12: #endif
14: /* data for runaway electron model */
15: typedef struct REctx_struct {
16: PetscErrorCode (*test)(TS, Vec, PetscInt, PetscReal, PetscBool, LandauCtx *, struct REctx_struct *);
17: PetscErrorCode (*impuritySrcRate)(PetscReal, PetscReal *, LandauCtx *);
18: PetscErrorCode (*E)(Vec, Vec, PetscInt, PetscReal, LandauCtx *, PetscReal *);
19: PetscReal T_cold; /* temperature of newly ionized electrons and impurity ions */
20: PetscReal ion_potential; /* ionization potential of impurity */
21: PetscReal Ne_ion; /* effective number of electrons shed in ioization of impurity */
22: PetscReal Ez_initial;
23: PetscReal L; /* inductance */
24: Vec X_0;
25: PetscInt imp_idx; /* index for impurity ionizing sink */
26: PetscReal pulse_start;
27: PetscReal pulse_width;
28: PetscReal pulse_rate;
29: PetscReal current_rate;
30: PetscInt plotIdx;
31: PetscInt plotStep;
32: PetscInt idx; /* cache */
33: PetscReal j; /* cache */
34: PetscReal plotDt;
35: PetscBool plotting;
36: PetscBool use_spitzer_eta;
37: PetscInt print_period;
38: PetscInt grid_view_idx;
39: } REctx;
41: static const PetscReal kev_joul = 6.241506479963235e+15; /* 1/1000e */
43: #define RE_CUT 3.
44: /* < v, u_re * v * q > */
45: static void f0_j_re(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
46: {
47: PetscReal n_e = PetscRealPart(u[0]);
48: if (dim == 2) {
49: if (x[1] > RE_CUT || x[1] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
50: *f0 = n_e * 2. * PETSC_PI * x[0] * x[1] * constants[0]; /* n * r * v_|| * q */
51: } else {
52: *f0 = 0;
53: }
54: } else {
55: if (x[2] > RE_CUT || x[2] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
56: *f0 = n_e * x[2] * constants[0];
57: } else {
58: *f0 = 0;
59: }
60: }
61: }
63: /* sum < v, u*v*q > */
64: static void f0_jz_sum(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar q[], PetscScalar *f0)
65: {
66: PetscInt ii;
67: f0[0] = 0;
68: if (dim == 2) {
69: for (ii = 0; ii < Nf; ii++) f0[0] += u[ii] * 2. * PETSC_PI * x[0] * x[1] * q[ii]; /* n * r * v_|| * q * v_0 */
70: } else {
71: for (ii = 0; ii < Nf; ii++) f0[0] += u[ii] * x[2] * q[ii]; /* n * v_|| * q * v_0 */
72: }
73: }
75: /* < v, n_e > */
76: static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
77: {
78: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
79: if (dim == 2) f0[0] = 2. * PETSC_PI * x[0] * u[ii];
80: else f0[0] = u[ii];
81: }
83: /* < v, n_e v_|| > */
84: static void f0_vz(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
85: {
86: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
87: if (dim == 2) f0[0] = u[ii] * 2. * PETSC_PI * x[0] * x[1]; /* n r v_|| */
88: else f0[0] = u[ii] * x[2]; /* n v_|| */
89: }
91: /* < v, n_e (v-shift) > */
92: static void f0_ve_shift(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
93: {
94: PetscReal vz = numConstants > 0 ? PetscRealPart(constants[0]) : 0;
95: if (dim == 2) *f0 = u[0] * 2. * PETSC_PI * x[0] * PetscSqrtReal(x[0] * x[0] + (x[1] - vz) * (x[1] - vz)); /* n r v */
96: else {
97: *f0 = u[0] * PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + (x[2] - vz) * (x[2] - vz)); /* n v */
98: }
99: }
101: /* CalculateE - Calculate the electric field */
102: /* T -- Electron temperature */
103: /* n -- Electron density */
104: /* lnLambda -- */
105: /* eps0 -- */
106: /* E -- output E, input \hat E */
107: static PetscReal CalculateE(PetscReal Tev, PetscReal n, PetscReal lnLambda, PetscReal eps0, PetscReal *E)
108: {
109: PetscReal c, e, m;
111: PetscFunctionBegin;
112: c = 299792458.0;
113: e = 1.602176e-19;
114: m = 9.10938e-31;
115: if (1) {
116: double Ec, Ehat = *E, betath = PetscSqrtReal(2 * Tev * e / (m * c * c)), j0 = Ehat * 7 / (PetscSqrtReal(2) * 2) * PetscPowReal(betath, 3) * n * e * c;
117: Ec = n * lnLambda * PetscPowReal(e, 3) / (4 * PETSC_PI * PetscPowReal(eps0, 2) * m * c * c);
118: *E = Ec;
119: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CalculateE j0=%g Ec = %g\n", j0, Ec));
120: } else {
121: PetscReal Ed, vth;
122: vth = PetscSqrtReal(8 * Tev * e / (m * PETSC_PI));
123: Ed = n * lnLambda * PetscPowReal(e, 3) / (4 * PETSC_PI * PetscPowReal(eps0, 2) * m * vth * vth);
124: *E = Ed;
125: }
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscReal Spitzer(PetscReal m_e, PetscReal e, PetscReal Z, PetscReal epsilon0, PetscReal lnLam, PetscReal kTe_joules)
130: {
131: PetscReal Fz = (1 + 1.198 * Z + 0.222 * Z * Z) / (1 + 2.966 * Z + 0.753 * Z * Z), eta;
132: eta = Fz * 4. / 3. * PetscSqrtReal(2. * PETSC_PI) * Z * PetscSqrtReal(m_e) * PetscSqr(e) * lnLam * PetscPowReal(4 * PETSC_PI * epsilon0, -2.) * PetscPowReal(kTe_joules, -1.5);
133: return eta;
134: }
136: static PetscErrorCode testNone(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
137: {
138: PetscFunctionBeginUser;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: static PetscErrorCode testSpitzer(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
143: {
144: PetscInt ii, nDMs;
145: PetscDS prob;
146: static PetscReal old_ratio = 1e10;
147: TSConvergedReason reason;
148: PetscReal J, J_re, spit_eta, Te_kev = 0, E, ratio, Z, n_e, v, v2;
149: PetscScalar user[2] = {0., ctx->charges[0]}, q[LANDAU_MAX_SPECIES], tt[LANDAU_MAX_SPECIES], vz;
150: PetscReal dt;
151: DM pack, plexe = ctx->plex[0], plexi = (ctx->num_grids == 1) ? NULL : ctx->plex[1];
152: Vec *XsubArray;
154: PetscFunctionBeginUser;
155: PetscCheck(ctx->num_species == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ctx->num_species %" PetscInt_FMT " != 2", ctx->num_species);
156: PetscCall(VecGetDM(X, &pack));
157: PetscCheck(pack, PETSC_COMM_SELF, PETSC_ERR_PLIB, "no DM");
158: PetscCall(DMCompositeGetNumberDM(pack, &nDMs));
159: PetscCheck(nDMs == ctx->num_grids * ctx->batch_sz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nDMs != ctx->num_grids*ctx->batch_sz %" PetscInt_FMT " != %" PetscInt_FMT, nDMs, ctx->num_grids * ctx->batch_sz);
160: PetscCall(PetscMalloc(sizeof(*XsubArray) * nDMs, &XsubArray));
161: PetscCall(DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray)); // read only
162: PetscCall(TSGetTimeStep(ts, &dt));
163: /* get current for each grid */
164: for (ii = 0; ii < ctx->num_species; ii++) q[ii] = ctx->charges[ii];
165: PetscCall(DMGetDS(plexe, &prob));
166: PetscCall(PetscDSSetConstants(prob, 2, &q[0]));
167: PetscCall(PetscDSSetObjective(prob, 0, &f0_jz_sum));
168: PetscCall(DMPlexComputeIntegralFEM(plexe, XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, 0)], tt, NULL));
169: J = -ctx->n_0 * ctx->v_0 * PetscRealPart(tt[0]);
170: if (plexi) { // add first (only) ion
171: PetscCall(DMGetDS(plexi, &prob));
172: PetscCall(PetscDSSetConstants(prob, 1, &q[1]));
173: PetscCall(PetscDSSetObjective(prob, 0, &f0_jz_sum));
174: PetscCall(DMPlexComputeIntegralFEM(plexi, XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, 1)], tt, NULL));
175: J += -ctx->n_0 * ctx->v_0 * PetscRealPart(tt[0]);
176: }
177: /* get N_e */
178: PetscCall(DMGetDS(plexe, &prob));
179: PetscCall(PetscDSSetConstants(prob, 1, user));
180: PetscCall(PetscDSSetObjective(prob, 0, &f0_n));
181: PetscCall(DMPlexComputeIntegralFEM(plexe, XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, 0)], tt, NULL));
182: n_e = PetscRealPart(tt[0]) * ctx->n_0;
183: /* Z */
184: Z = -ctx->charges[1] / ctx->charges[0];
185: /* remove drift */
186: if (0) {
187: user[0] = 0; // electrons
188: PetscCall(DMGetDS(plexe, &prob));
189: PetscCall(PetscDSSetConstants(prob, 1, user));
190: PetscCall(PetscDSSetObjective(prob, 0, &f0_vz));
191: PetscCall(DMPlexComputeIntegralFEM(plexe, XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, 0)], tt, NULL));
192: vz = ctx->n_0 * PetscRealPart(tt[0]) / n_e; /* non-dimensional */
193: } else vz = 0;
194: /* thermal velocity */
195: PetscCall(DMGetDS(plexe, &prob));
196: PetscCall(PetscDSSetConstants(prob, 1, &vz));
197: PetscCall(PetscDSSetObjective(prob, 0, &f0_ve_shift));
198: PetscCall(DMPlexComputeIntegralFEM(plexe, XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, 0)], tt, NULL));
199: v = ctx->n_0 * ctx->v_0 * PetscRealPart(tt[0]) / n_e; /* remove number density to get velocity */
200: v2 = PetscSqr(v); /* use real space: m^2 / s^2 */
201: Te_kev = (v2 * ctx->masses[0] * PETSC_PI / 8) * kev_joul; /* temperature in kev */
202: spit_eta = Spitzer(ctx->masses[0], -ctx->charges[0], Z, ctx->epsilon0, ctx->lambdas[0][1], Te_kev / kev_joul); /* kev --> J (kT) */
203: if (0) {
204: PetscCall(DMGetDS(plexe, &prob));
205: PetscCall(PetscDSSetConstants(prob, 1, q));
206: PetscCall(PetscDSSetObjective(prob, 0, &f0_j_re));
207: PetscCall(DMPlexComputeIntegralFEM(plexe, XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, 0)], tt, NULL));
208: } else tt[0] = 0;
209: J_re = -ctx->n_0 * ctx->v_0 * PetscRealPart(tt[0]);
210: PetscCall(DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray)); // read only
211: PetscCall(PetscFree(XsubArray));
213: if (rectx->use_spitzer_eta) {
214: E = ctx->Ez = spit_eta * (rectx->j - J_re);
215: } else {
216: E = ctx->Ez; /* keep real E */
217: rectx->j = J; /* cache */
218: }
220: ratio = E / J / spit_eta;
221: if (stepi > 10 && !rectx->use_spitzer_eta && (old_ratio - ratio < 1.e-6)) {
222: rectx->pulse_start = time + 0.98 * dt;
223: rectx->use_spitzer_eta = PETSC_TRUE;
224: }
225: PetscCall(TSGetConvergedReason(ts, &reason));
226: PetscCall(TSGetConvergedReason(ts, &reason));
227: if (rectx->plotting || stepi == 0 || reason || rectx->pulse_start == time + 0.98 * dt) {
228: PetscCall(PetscPrintf(ctx->comm, "testSpitzer: %4" PetscInt_FMT ") time=%11.4e n_e= %10.3e E= %10.3e J= %10.3e J_re= %10.3e %.3g%% Te_kev= %10.3e Z_eff=%g E/J to eta ratio= %g (diff=%g) %s %s spit_eta=%g\n", stepi, (double)time,
229: (double)(n_e / ctx->n_0), (double)ctx->Ez, (double)J, (double)J_re, (double)(100 * J_re / J), (double)Te_kev, (double)Z, (double)ratio, (double)(old_ratio - ratio), rectx->use_spitzer_eta ? "using Spitzer eta*J E" : "constant E", rectx->pulse_start != time + 0.98 * dt ? "normal" : "transition", (double)spit_eta));
230: PetscCheck(rectx->pulse_start != (time + 0.98 * dt), PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Spitzer complete ratio=%g", (double)ratio);
231: }
232: old_ratio = ratio;
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: static const double ppp = 2;
237: static void f0_0_diff_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
238: {
239: LandauCtx *ctx = (LandauCtx *)constants;
240: REctx *rectx = (REctx *)ctx->data;
241: PetscInt ii = rectx->idx, i;
242: const PetscReal kT_m = ctx->k * ctx->thermal_temps[ii] / ctx->masses[ii]; /* kT/m */
243: const PetscReal n = ctx->n[ii];
244: PetscReal diff, f_maxwell, v2 = 0, theta = 2 * kT_m / (ctx->v_0 * ctx->v_0); /* theta = 2kT/mc^2 */
245: for (i = 0; i < dim; ++i) v2 += x[i] * x[i];
246: f_maxwell = n * PetscPowReal(PETSC_PI * theta, -1.5) * (PetscExpReal(-v2 / theta));
247: diff = 2. * PETSC_PI * x[0] * (PetscRealPart(u[ii]) - f_maxwell);
248: f0[0] = PetscPowReal(diff, ppp);
249: }
250: static void f0_0_maxwellian_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
251: {
252: LandauCtx *ctx = (LandauCtx *)constants;
253: REctx *rectx = (REctx *)ctx->data;
254: PetscInt ii = rectx->idx, i;
255: const PetscReal kT_m = ctx->k * ctx->thermal_temps[ii] / ctx->masses[ii]; /* kT/m */
256: const PetscReal n = ctx->n[ii];
257: PetscReal f_maxwell, v2 = 0, theta = 2 * kT_m / (ctx->v_0 * ctx->v_0); /* theta = 2kT/mc^2 */
258: for (i = 0; i < dim; ++i) v2 += x[i] * x[i];
259: f_maxwell = 2. * PETSC_PI * x[0] * n * PetscPowReal(PETSC_PI * theta, -1.5) * (PetscExpReal(-v2 / theta));
260: f0[0] = PetscPowReal(f_maxwell, ppp);
261: }
263: static PetscErrorCode testStable(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
264: {
265: PetscDS prob;
266: Vec X2;
267: PetscReal ediff, idiff = 0, lpm0, lpm1 = 1;
268: PetscScalar tt[LANDAU_MAX_SPECIES];
269: DM dm, plex = ctx->plex[0];
271: PetscFunctionBeginUser;
272: PetscCall(VecGetDM(X, &dm));
273: PetscCall(DMGetDS(plex, &prob));
274: PetscCall(VecDuplicate(X, &X2));
275: PetscCall(VecCopy(X, X2));
276: if (!rectx->X_0) {
277: PetscCall(VecDuplicate(X, &rectx->X_0));
278: PetscCall(VecCopy(X, rectx->X_0));
279: }
280: PetscCall(VecAXPY(X, -1.0, rectx->X_0));
281: PetscCall(PetscDSSetConstants(prob, sizeof(LandauCtx) / sizeof(PetscScalar), (PetscScalar *)ctx));
282: rectx->idx = 0;
283: PetscCall(PetscDSSetObjective(prob, 0, &f0_0_diff_lp));
284: PetscCall(DMPlexComputeIntegralFEM(plex, X2, tt, NULL));
285: ediff = PetscPowReal(PetscRealPart(tt[0]), 1. / ppp);
286: PetscCall(PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp));
287: PetscCall(DMPlexComputeIntegralFEM(plex, X2, tt, NULL));
288: lpm0 = PetscPowReal(PetscRealPart(tt[0]), 1. / ppp);
289: if (ctx->num_species > 1) {
290: rectx->idx = 1;
291: PetscCall(PetscDSSetObjective(prob, 0, &f0_0_diff_lp));
292: PetscCall(DMPlexComputeIntegralFEM(plex, X2, tt, NULL));
293: idiff = PetscPowReal(PetscRealPart(tt[0]), 1. / ppp);
294: PetscCall(PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp));
295: PetscCall(DMPlexComputeIntegralFEM(plex, X2, tt, NULL));
296: lpm1 = PetscPowReal(PetscRealPart(tt[0]), 1. / ppp);
297: }
298: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s %" PetscInt_FMT ") time=%10.3e n-%d norm electrons/max=%20.13e ions/max=%20.13e\n", "----", stepi, (double)time, (int)ppp, (double)(ediff / lpm0), (double)(idiff / lpm1)));
299: /* view */
300: PetscCall(VecCopy(X2, X));
301: PetscCall(VecDestroy(&X2));
302: if (islast) {
303: PetscCall(VecDestroy(&rectx->X_0));
304: rectx->X_0 = NULL;
305: }
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode EInduction(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
310: {
311: REctx *rectx = (REctx *)ctx->data;
312: PetscInt ii;
313: DM dm, plex;
314: PetscScalar tt[LANDAU_MAX_SPECIES], qv0[LANDAU_MAX_SPECIES];
315: PetscReal dJ_dt;
316: PetscDS prob;
318: PetscFunctionBeginUser;
319: for (ii = 0; ii < ctx->num_species; ii++) qv0[ii] = ctx->charges[ii] * ctx->v_0;
320: PetscCall(VecGetDM(X, &dm));
321: PetscCall(DMGetDS(dm, &prob));
322: PetscCall(DMConvert(dm, DMPLEX, &plex));
323: /* get d current / dt */
324: PetscCall(PetscDSSetConstants(prob, ctx->num_species, qv0));
325: PetscCall(PetscDSSetObjective(prob, 0, &f0_jz_sum));
326: PetscCheck(X_t, PETSC_COMM_SELF, PETSC_ERR_PLIB, "X_t");
327: PetscCall(DMPlexComputeIntegralFEM(plex, X_t, tt, NULL));
328: dJ_dt = -ctx->n_0 * PetscRealPart(tt[0]) / ctx->t_0;
329: /* E induction */
330: *a_E = -rectx->L * dJ_dt + rectx->Ez_initial;
331: PetscCall(DMDestroy(&plex));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: static PetscErrorCode EConstant(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
336: {
337: PetscFunctionBeginUser;
338: *a_E = ctx->Ez;
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode ENone(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
343: {
344: PetscFunctionBeginUser;
345: *a_E = 0;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /* ------------------------------------------------------------------- */
350: /*
351: FormSource - Evaluates source terms F(t).
353: Input Parameters:
354: . ts - the TS context
355: . time -
356: . X_dummmy - input vector
357: . dummy - optional user-defined context, as set by SNESSetFunction()
359: Output Parameter:
360: . F - function vector
361: */
362: static PetscErrorCode FormSource(TS ts, PetscReal ftime, Vec X_dummmy, Vec F, void *dummy)
363: {
364: PetscReal new_imp_rate;
365: LandauCtx *ctx;
366: DM pack;
367: REctx *rectx;
369: PetscFunctionBeginUser;
370: PetscCall(TSGetDM(ts, &pack));
371: PetscCall(DMGetApplicationContext(pack, &ctx));
372: rectx = (REctx *)ctx->data;
373: /* check for impurities */
374: PetscCall(rectx->impuritySrcRate(ftime, &new_imp_rate, ctx));
375: if (new_imp_rate != 0) {
376: if (new_imp_rate != rectx->current_rate) {
377: PetscInt ii;
378: PetscReal dne_dt, dni_dt, tilda_ns[LANDAU_MAX_SPECIES], temps[LANDAU_MAX_SPECIES];
379: Vec globFarray[LANDAU_MAX_GRIDS * LANDAU_MAX_BATCH_SZ];
380: rectx->current_rate = new_imp_rate;
381: for (ii = 1; ii < LANDAU_MAX_SPECIES; ii++) tilda_ns[ii] = 0;
382: for (ii = 1; ii < LANDAU_MAX_SPECIES; ii++) temps[ii] = 1;
383: dni_dt = new_imp_rate /* *ctx->t_0 */; /* fully ionized immediately, no normalize, stay in non-dim */
384: dne_dt = new_imp_rate * rectx->Ne_ion /* *ctx->t_0 */;
385: tilda_ns[0] = dne_dt;
386: tilda_ns[rectx->imp_idx] = dni_dt;
387: temps[0] = rectx->T_cold;
388: temps[rectx->imp_idx] = rectx->T_cold;
389: PetscCall(PetscInfo(ctx->plex[0], "\tHave new_imp_rate= %10.3e time= %10.3e de/dt= %10.3e di/dt= %10.3e ***\n", (double)new_imp_rate, (double)ftime, (double)dne_dt, (double)dni_dt));
390: PetscCall(DMCompositeGetAccessArray(pack, F, ctx->num_grids * ctx->batch_sz, NULL, globFarray));
391: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
392: /* add it */
393: PetscCall(DMPlexLandauAddMaxwellians(ctx->plex[grid], globFarray[LAND_PACK_IDX(0, grid)], ftime, temps, tilda_ns, grid, 0, 1, ctx));
394: }
395: // Does DMCompositeRestoreAccessArray copy the data back? (no)
396: PetscCall(DMCompositeRestoreAccessArray(pack, F, ctx->num_grids * ctx->batch_sz, NULL, globFarray));
397: }
398: } else {
399: PetscCall(VecZeroEntries(F));
400: rectx->current_rate = 0;
401: }
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx)
406: {
407: LandauCtx *ctx = (LandauCtx *)actx; /* user-defined application context */
408: REctx *rectx = (REctx *)ctx->data;
409: DM pack = NULL;
410: Vec globXArray[LANDAU_MAX_GRIDS * LANDAU_MAX_BATCH_SZ];
411: TSConvergedReason reason;
413: PetscFunctionBeginUser;
414: PetscCall(TSGetConvergedReason(ts, &reason));
415: if (rectx->grid_view_idx != -1 || (reason && ctx->verbose > 3)) {
416: PetscCall(VecGetDM(X, &pack));
417: PetscCall(DMCompositeGetAccessArray(pack, X, ctx->num_grids * ctx->batch_sz, NULL, globXArray));
418: }
419: if (stepi > rectx->plotStep && rectx->plotting) {
420: rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */
421: rectx->plotIdx++;
422: }
423: /* view */
424: if (time / rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) {
425: if ((reason || stepi == 0 || rectx->plotIdx % rectx->print_period == 0) && ctx->verbose > 1) {
426: /* print norms */
427: PetscCall(DMPlexLandauPrintNorms(X, stepi));
428: }
429: if (!rectx->plotting) { /* first step of possible backtracks */
430: rectx->plotting = PETSC_TRUE;
431: /* diagnostics + change E field with Sptizer (not just a monitor) */
432: PetscCall(rectx->test(ts, X, stepi, time, reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx));
433: } else {
434: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t\t ERROR SKIP test spit ------\n"));
435: rectx->plotting = PETSC_TRUE;
436: }
437: if (rectx->grid_view_idx != -1) {
438: PetscCall(PetscObjectSetName((PetscObject)globXArray[LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx)], rectx->grid_view_idx == 0 ? "ue" : "ui"));
439: /* view, overwrite step when back tracked */
440: PetscCall(DMSetOutputSequenceNumber(ctx->plex[rectx->grid_view_idx], rectx->plotIdx, time * ctx->t_0));
441: PetscCall(VecViewFromOptions(globXArray[LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx)], NULL, "-ex2_vec_view"));
442: }
443: rectx->plotStep = stepi;
444: } else {
445: if (rectx->plotting) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ERROR rectx->plotting=%s step %" PetscInt_FMT "\n", PetscBools[rectx->plotting], stepi));
446: /* diagnostics + change E field with Sptizer (not just a monitor) - can we lag this? */
447: PetscCall(rectx->test(ts, X, stepi, time, reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx));
448: }
449: /* parallel check that only works of all batches are identical */
450: if (reason && ctx->verbose > 3 && ctx->batch_sz > 1) {
451: PetscReal val, rval;
452: PetscMPIInt rank;
453: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
454: for (PetscInt grid = 0; grid < ctx->num_grids; grid++) {
455: PetscInt nerrors = 0;
456: for (PetscInt i = 0; i < ctx->batch_sz; i++) {
457: PetscCall(VecNorm(globXArray[LAND_PACK_IDX(i, grid)], NORM_2, &val));
458: if (i == 0) rval = val;
459: else if ((val = PetscAbs(val - rval) / rval) > 1000 * PETSC_MACHINE_EPSILON) {
460: PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] Warning %" PetscInt_FMT ".%" PetscInt_FMT ") diff = %2.15e\n", rank, grid, i, (double)val));
461: nerrors++;
462: }
463: }
464: if (nerrors) {
465: PetscCall(PetscPrintf(PETSC_COMM_SELF, " ***** [%d] ERROR max %" PetscInt_FMT " errors\n", rank, nerrors));
466: } else {
467: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%d] %" PetscInt_FMT ") batch consistency check OK\n", rank, grid));
468: }
469: }
470: }
471: rectx->idx = 0;
472: if (rectx->grid_view_idx != -1 || (reason && ctx->verbose > 3)) PetscCall(DMCompositeRestoreAccessArray(pack, X, ctx->num_grids * ctx->batch_sz, NULL, globXArray));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: PetscErrorCode PreStep(TS ts)
477: {
478: LandauCtx *ctx;
479: REctx *rectx;
480: DM dm;
481: PetscInt stepi;
482: PetscReal time;
483: Vec X;
485: PetscFunctionBeginUser;
486: /* not used */
487: PetscCall(TSGetDM(ts, &dm));
488: PetscCall(TSGetTime(ts, &time));
489: PetscCall(TSGetSolution(ts, &X));
490: PetscCall(DMGetApplicationContext(dm, &ctx));
491: rectx = (REctx *)ctx->data;
492: PetscCall(TSGetStepNumber(ts, &stepi));
493: /* update E */
494: PetscCall(rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez));
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /* model for source of non-ionized impurities, profile provided by model, in du/dt form in normalized units (tricky because n_0 is normalized with electrons) */
499: static PetscErrorCode stepSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
500: {
501: REctx *rectx = (REctx *)ctx->data;
503: PetscFunctionBeginUser;
504: if (time >= rectx->pulse_start) *rho = rectx->pulse_rate;
505: else *rho = 0.;
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
508: static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
509: {
510: PetscFunctionBeginUser;
511: *rho = 0.;
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
514: static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
515: {
516: REctx *rectx = (REctx *)ctx->data;
518: PetscFunctionBeginUser;
519: PetscCheck(rectx->pulse_start != PETSC_MAX_REAL, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "'-ex2_pulse_start_time X' must be used with '-ex2_impurity_source_type pulse'");
520: if (time < rectx->pulse_start || time > rectx->pulse_start + 3 * rectx->pulse_width) *rho = 0;
521: else {
522: double x = PetscSinReal((time - rectx->pulse_start) / (3 * rectx->pulse_width) * 2 * PETSC_PI - PETSC_PI / 2) + 1; /* 0:2, integrates to 1.0 */
523: *rho = rectx->pulse_rate * x / (3 * rectx->pulse_width);
524: if (!rectx->use_spitzer_eta) rectx->use_spitzer_eta = PETSC_TRUE; /* use it next time */
525: }
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
531: static PetscErrorCode ProcessREOptions(REctx *rectx, const LandauCtx *ctx, DM dm, const char prefix[])
532: {
533: PetscFunctionList plist = NULL, testlist = NULL, elist = NULL;
534: char pname[256], testname[256], ename[256];
535: DM dm_dummy;
536: PetscBool Connor_E = PETSC_FALSE;
538: PetscFunctionBeginUser;
539: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm_dummy));
540: rectx->Ne_ion = 1; /* number of electrons given up by impurity ion */
541: rectx->T_cold = .005; /* kev */
542: rectx->ion_potential = 15; /* ev */
543: rectx->L = 2;
544: rectx->X_0 = NULL;
545: rectx->imp_idx = ctx->num_species - 1; /* default ionized impurity as last one */
546: rectx->pulse_start = PETSC_MAX_REAL;
547: rectx->pulse_width = 1;
548: rectx->plotStep = PETSC_INT_MAX;
549: rectx->pulse_rate = 1.e-1;
550: rectx->current_rate = 0;
551: rectx->plotIdx = 0;
552: rectx->j = 0;
553: rectx->plotDt = 1.0;
554: rectx->plotting = PETSC_FALSE;
555: rectx->use_spitzer_eta = PETSC_FALSE;
556: rectx->idx = 0;
557: rectx->print_period = 10;
558: rectx->grid_view_idx = -1; // do not get if not needed
559: /* Register the available impurity sources */
560: PetscCall(PetscFunctionListAdd(&plist, "step", &stepSrc));
561: PetscCall(PetscFunctionListAdd(&plist, "none", &zeroSrc));
562: PetscCall(PetscFunctionListAdd(&plist, "pulse", &pulseSrc));
563: PetscCall(PetscStrncpy(pname, "none", sizeof(pname)));
564: PetscCall(PetscFunctionListAdd(&testlist, "none", &testNone));
565: PetscCall(PetscFunctionListAdd(&testlist, "spitzer", &testSpitzer));
566: PetscCall(PetscFunctionListAdd(&testlist, "stable", &testStable));
567: PetscCall(PetscStrncpy(testname, "none", sizeof(testname)));
568: PetscCall(PetscFunctionListAdd(&elist, "none", &ENone));
569: PetscCall(PetscFunctionListAdd(&elist, "induction", &EInduction));
570: PetscCall(PetscFunctionListAdd(&elist, "constant", &EConstant));
571: PetscCall(PetscStrncpy(ename, "constant", sizeof(ename)));
573: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");
574: PetscCall(PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "ex2.c", rectx->plotDt, &rectx->plotDt, NULL));
575: if (rectx->plotDt < 0) rectx->plotDt = 1e30;
576: if (rectx->plotDt == 0) rectx->plotDt = 1e-30;
577: PetscCall(PetscOptionsInt("-ex2_print_period", "Plotting interval", "ex2.c", rectx->print_period, &rectx->print_period, NULL));
578: PetscCall(PetscOptionsInt("-ex2_grid_view_idx", "grid_view_idx", "ex2.c", rectx->grid_view_idx, &rectx->grid_view_idx, NULL));
579: PetscCheck(rectx->grid_view_idx < ctx->num_grids || rectx->grid_view_idx == -1, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "rectx->grid_view_idx (%" PetscInt_FMT ") >= ctx->num_grids (%" PetscInt_FMT ")", rectx->imp_idx, ctx->num_grids);
580: PetscCall(PetscOptionsFList("-ex2_impurity_source_type", "Name of impurity source to run", "", plist, pname, pname, sizeof(pname), NULL));
581: PetscCall(PetscOptionsFList("-ex2_test_type", "Name of test to run", "", testlist, testname, testname, sizeof(testname), NULL));
582: PetscCall(PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL));
583: PetscCheck((rectx->imp_idx < ctx->num_species && rectx->imp_idx >= 1) || ctx->num_species <= 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "index of sink for impurities ions is out of range (%" PetscInt_FMT "), must be > 0 && < NS", rectx->imp_idx);
584: PetscCall(PetscOptionsFList("-ex2_e_field_type", "Electric field type", "", elist, ename, ename, sizeof(ename), NULL));
585: rectx->Ne_ion = -ctx->charges[rectx->imp_idx] / ctx->charges[0];
586: PetscCall(PetscOptionsReal("-ex2_t_cold", "Temperature of cold electron and ions after ionization in keV", "none", rectx->T_cold, &rectx->T_cold, NULL));
587: PetscCall(PetscOptionsReal("-ex2_pulse_start_time", "Time at which pulse happens for 'pulse' source", "none", rectx->pulse_start, &rectx->pulse_start, NULL));
588: PetscCall(PetscOptionsReal("-ex2_pulse_width_time", "Width of pulse 'pulse' source", "none", rectx->pulse_width, &rectx->pulse_width, NULL));
589: PetscCall(PetscOptionsReal("-ex2_pulse_rate", "Number density of pulse for 'pulse' source", "none", rectx->pulse_rate, &rectx->pulse_rate, NULL));
590: rectx->T_cold *= 1.16e7; /* convert to Kelvin */
591: PetscCall(PetscOptionsReal("-ex2_ion_potential", "Potential to ionize impurity (should be array) in ev", "none", rectx->ion_potential, &rectx->ion_potential, NULL));
592: PetscCall(PetscOptionsReal("-ex2_inductance", "Inductance E field", "none", rectx->L, &rectx->L, NULL));
593: PetscCall(PetscOptionsBool("-ex2_connor_e_field_units", "Scale Ex but Connor-Hastie E_c", "none", Connor_E, &Connor_E, NULL));
594: PetscCall(PetscInfo(dm_dummy, "Num electrons from ions=%g, T_cold=%10.3e, ion potential=%10.3e, E_z=%10.3e v_0=%10.3e\n", (double)rectx->Ne_ion, (double)rectx->T_cold, (double)rectx->ion_potential, (double)ctx->Ez, (double)ctx->v_0));
595: PetscOptionsEnd();
596: /* get impurity source rate function */
597: PetscCall(PetscFunctionListFind(plist, pname, &rectx->impuritySrcRate));
598: PetscCheck(rectx->impuritySrcRate, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "No impurity source function found '%s'", pname);
599: PetscCall(PetscFunctionListFind(testlist, testname, &rectx->test));
600: PetscCheck(rectx->test, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "No test found '%s'", testname);
601: PetscCall(PetscFunctionListFind(elist, ename, &rectx->E));
602: PetscCheck(rectx->E, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "No E field function found '%s'", ename);
603: PetscCall(PetscFunctionListDestroy(&plist));
604: PetscCall(PetscFunctionListDestroy(&testlist));
605: PetscCall(PetscFunctionListDestroy(&elist));
607: /* convert E from Connor-Hastie E_c units to real if doing Spitzer E */
608: if (Connor_E) {
609: PetscReal E = ctx->Ez, Tev = ctx->thermal_temps[0] * 8.621738e-5, n = ctx->n_0 * ctx->n[0];
610: CalculateE(Tev, n, ctx->lambdas[0][1], ctx->epsilon0, &E);
611: ((LandauCtx *)ctx)->Ez *= E;
612: }
613: PetscCall(DMDestroy(&dm_dummy));
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: int main(int argc, char **argv)
618: {
619: DM pack;
620: Vec X;
621: PetscInt dim = 2, nDMs;
622: TS ts;
623: Mat J;
624: PetscDS prob;
625: LandauCtx *ctx;
626: REctx *rectx;
627: PetscMPIInt rank;
628: PetscLogStage stage;
630: PetscFunctionBeginUser;
631: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
632: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
633: if (rank) { /* turn off output stuff for duplicate runs */
634: PetscCall(PetscOptionsClearValue(NULL, "-ex2_dm_view"));
635: PetscCall(PetscOptionsClearValue(NULL, "-ex2_vec_view"));
636: PetscCall(PetscOptionsClearValue(NULL, "-ex2_vec_view_init"));
637: PetscCall(PetscOptionsClearValue(NULL, "-ex2_dm_view_init"));
638: PetscCall(PetscOptionsClearValue(NULL, "-info")); /* this does not work */
639: }
640: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
641: /* Create a mesh */
642: PetscCall(DMPlexLandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, "", &X, &J, &pack));
643: PetscCall(DMCompositeGetNumberDM(pack, &nDMs));
644: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
645: PetscCall(PetscObjectSetName((PetscObject)X, "f"));
646: PetscCall(DMGetApplicationContext(pack, &ctx));
647: PetscCall(DMSetUp(pack));
648: /* context */
649: PetscCall(PetscNew(&rectx));
650: ctx->data = rectx;
651: PetscCall(ProcessREOptions(rectx, ctx, pack, ""));
652: PetscCall(DMGetDS(pack, &prob));
653: if (rectx->grid_view_idx != -1) {
654: Vec *XsubArray = NULL;
655: PetscCall(PetscMalloc(sizeof(*XsubArray) * nDMs, &XsubArray));
656: PetscCall(DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray)); // read only
657: PetscCall(PetscObjectSetName((PetscObject)XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx)], rectx->grid_view_idx == 0 ? "ue" : "ui"));
658: PetscCall(DMSetOutputSequenceNumber(ctx->plex[rectx->grid_view_idx], 0, 0.0));
659: PetscCall(DMViewFromOptions(ctx->plex[rectx->grid_view_idx], NULL, "-ex2_dm_view"));
660: PetscCall(DMViewFromOptions(ctx->plex[rectx->grid_view_idx], NULL, "-ex2_dm_view_init"));
661: PetscCall(VecViewFromOptions(XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx)], NULL, "-ex2_vec_view")); // initial condition (monitor plots after step)
662: PetscCall(VecViewFromOptions(XsubArray[LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx)], NULL, "-ex2_vec_view_init")); // initial condition (monitor plots after step)
663: PetscCall(DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray)); // read only
664: PetscCall(PetscFree(XsubArray));
665: }
666: /* Create timestepping solver context */
667: PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
668: PetscCall(TSSetDM(ts, pack));
669: PetscCall(TSSetIFunction(ts, NULL, DMPlexLandauIFunction, NULL));
670: PetscCall(TSSetIJacobian(ts, J, J, DMPlexLandauIJacobian, NULL));
671: PetscCall(TSSetRHSFunction(ts, NULL, FormSource, NULL));
672: PetscCall(TSSetFromOptions(ts));
673: PetscCall(TSSetSolution(ts, X));
674: PetscCall(TSSetApplicationContext(ts, ctx));
675: PetscCall(TSMonitorSet(ts, Monitor, ctx, NULL));
676: PetscCall(TSSetPreStep(ts, PreStep));
677: rectx->Ez_initial = ctx->Ez; /* cache for induction calculation - applied E field */
678: if (1) { /* warm up an test just DMPlexLandauIJacobian */
679: Vec vec;
680: PetscInt nsteps;
681: PetscReal dt;
682: PetscCall(PetscLogStageRegister("Warmup", &stage));
683: PetscCall(PetscLogStagePush(stage));
684: PetscCall(VecDuplicate(X, &vec));
685: PetscCall(VecCopy(X, vec));
686: PetscCall(TSGetMaxSteps(ts, &nsteps));
687: PetscCall(TSGetTimeStep(ts, &dt));
688: PetscCall(TSSetMaxSteps(ts, 1));
689: PetscCall(TSSolve(ts, X));
690: PetscCall(TSSetMaxSteps(ts, nsteps));
691: PetscCall(TSSetStepNumber(ts, 0));
692: PetscCall(TSSetTime(ts, 0));
693: PetscCall(TSSetTimeStep(ts, dt));
694: rectx->plotIdx = 0;
695: rectx->plotting = PETSC_FALSE;
696: PetscCall(PetscLogStagePop());
697: PetscCall(VecCopy(vec, X));
698: PetscCall(VecDestroy(&vec));
699: PetscCall(PetscObjectStateIncrease((PetscObject)ctx->J));
700: }
701: /* go */
702: PetscCall(PetscLogStageRegister("Solve", &stage));
703: ctx->stage = 0; // lets not use this stage
704: PetscCall(PetscLogStagePush(stage));
705: #if defined(PETSC_HAVE_CUDA_NVTX)
706: nvtxRangePushA("ex2-TSSolve-warm");
707: #endif
708: PetscCall(TSSolve(ts, X));
709: #if defined(PETSC_HAVE_CUDA_NVTX)
710: nvtxRangePop();
711: #endif
712: PetscCall(PetscLogStagePop());
713: /* clean up */
714: PetscCall(DMPlexLandauDestroyVelocitySpace(&pack));
715: PetscCall(TSDestroy(&ts));
716: PetscCall(VecDestroy(&X));
717: PetscCall(PetscFree(rectx));
718: PetscCall(PetscFinalize());
719: return 0;
720: }
722: /*TEST
724: testset:
725: requires: p4est !complex double defined(PETSC_USE_DMLANDAU_2D)
726: output_file: output/ex2_0.out
727: args: -dm_landau_num_species_grid 1,1 -dm_landau_Ez 0 -petscspace_degree 3 -petscspace_poly_tensor 1 -dm_landau_type p4est -dm_landau_ion_masses 2 -dm_landau_ion_charges 1 -dm_landau_thermal_temps 5,5 -dm_landau_n 2,2 -dm_landau_n_0 5e19 -ts_monitor -snes_rtol 1.e-9 -snes_stol 1.e-14 -snes_monitor -snes_converged_reason -snes_max_it 10 -ts_type arkimex -ts_arkimex_type 1bee -ts_max_snes_failures unlimited -ts_rtol 1e-3 -ts_dt 1.e-2 -ts_max_time 1 -ts_adapt_clip .5,1.25 -ts_max_steps 2 -ts_adapt_scale_solve_failed 0.75 -ts_adapt_time_step_increase_delay 5 -dm_landau_amr_levels_max 2,2 -dm_landau_amr_re_levels 2 -dm_landau_re_radius 0 -ex2_impurity_source_type pulse -ex2_pulse_start_time 1e-1 -ex2_pulse_width_time 10 -ex2_pulse_rate 1e-2 -ex2_t_cold .05 -ex2_plot_dt 1e-1 -dm_refine 0 -dm_landau_gpu_assembly true -dm_landau_batch_size 2 -dm_landau_verbose 2 -dm_landau_domain_radius 5.,5.
728: test:
729: suffix: cpu
730: args: -dm_landau_device_type cpu -ksp_type bicg -pc_type jacobi
731: test:
732: suffix: kokkos
733: requires: kokkos_kernels !defined(PETSC_HAVE_CUDA_CLANG)
734: args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos -ksp_type bicg -pc_type jacobi
735: test:
736: suffix: kokkos_batch
737: requires: kokkos_kernels
738: args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos -ksp_type preonly -pc_type bjkokkos -pc_bjkokkos_ksp_type bicg -pc_bjkokkos_pc_type jacobi
739: test:
740: suffix: kokkos_batch_tfqmr
741: requires: kokkos_kernels !defined(PETSC_HAVE_CUDA_CLANG)
742: args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos -ksp_type preonly -pc_type bjkokkos -pc_bjkokkos_ksp_type tfqmr -pc_bjkokkos_pc_type jacobi
744: test:
745: requires: !complex double defined(PETSC_USE_DMLANDAU_2D) !cuda
746: suffix: single
747: nsize: 1
748: args: -dm_refine 2 -dm_landau_num_species_grid 1 -dm_landau_thermal_temps 1 -dm_landau_electron_shift 1.25 -petscspace_degree 3 -snes_converged_reason -ts_type beuler -ts_dt .1 -ex2_plot_dt .1 -ts_max_steps 1 -ex2_grid_view_idx 0 -ex2_dm_view -snes_rtol 1.e-13 -snes_stol 1.e-13 -dm_landau_verbose 2 -ex2_print_period 1 -ksp_type preonly -pc_type lu -dm_landau_device_type cpu -dm_landau_use_relativistic_corrections
750: testset:
751: requires: !complex double defined(PETSC_USE_DMLANDAU_2D)
752: nsize: 1
753: output_file: output/ex2_simplex.out
754: args: -dim 2 -dm_landau_num_species_grid 1,1 -petscspace_degree 2 -dm_landau_simplex -dm_landau_ion_masses 2 -dm_landau_ion_charges 1 -dm_landau_thermal_temps 2,1 -dm_landau_n 1,1 -snes_rtol 1e-15 -snes_stol 1e-15 -snes_monitor -ts_type beuler -snes_converged_reason -ts_exact_final_time stepover -ts_dt .1 -ts_max_steps 1 -ts_max_snes_failures unlimited -ksp_type preonly -pc_type lu -dm_landau_verbose 2 -ex2_grid_view_idx 0 -ex2_dm_view -dm_refine 1 -ksp_type bicg -pc_type jacobi
755: test:
756: suffix: simplex
757: args: -dm_landau_device_type cpu
758: test:
759: suffix: simplexkokkos
760: requires: kokkos_kernels !defined(PETSC_HAVE_CUDA_CLANG) !sycl
761: args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos
763: test:
764: requires: double !defined(PETSC_USE_DMLANDAU_2D)
765: suffix: sphere_3d
766: nsize: 1
767: args: -dim 3 -dm_landau_thermal_temps 2 -petscspace_degree 2 -ts_type beuler -ts_dt .1 -ts_max_steps 1 -dm_landau_verbose 2 -ksp_type preonly -pc_type lu -dm_landau_device_type cpu -snes_rtol 1.e-14 -snes_stol 1.e-14 -snes_converged_reason -dm_landau_sphere -dm_landau_sphere_inner_radius_90degree_scale .55 -dm_landau_sphere_inner_radius_45degree_scale .5
769: TEST*/