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