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>

  9: /* data for runaway electron model */
 10: typedef struct REctx_struct {
 11:   PetscErrorCode (*test)(TS, Vec, PetscInt, PetscReal, PetscBool,  LandauCtx *, struct REctx_struct *);
 12:   PetscErrorCode (*impuritySrcRate)(PetscReal, PetscReal *, LandauCtx*);
 13:   PetscErrorCode (*E)(Vec, Vec, PetscInt, PetscReal, LandauCtx*, PetscReal *);
 14:   PetscReal     T_cold;        /* temperature of newly ionized electrons and impurity ions */
 15:   PetscReal     ion_potential; /* ionization potential of impurity */
 16:   PetscReal     Ne_ion;        /* effective number of electrons shed in ioization of impurity */
 17:   PetscReal     Ez_initial;
 18:   PetscReal     L;             /* inductance */
 19:   Vec           X_0;
 20:   PetscInt      imp_idx;       /* index for impurity ionizing sink */
 21:   PetscReal     pulse_start;
 22:   PetscReal     pulse_width;
 23:   PetscReal     pulse_rate;
 24:   PetscReal     current_rate;
 25:   PetscInt      plotIdx;
 26:   PetscInt      plotStep;
 27:   PetscInt      idx; /* cache */
 28:   PetscReal     j; /* cache */
 29:   PetscReal     plotDt;
 30:   PetscBool     plotting;
 31:   PetscBool     use_spitzer_eta;
 32:   PetscInt      print_period;
 33:   PetscInt      grid_view_idx;
 34: } REctx;

 36: static const PetscReal kev_joul = 6.241506479963235e+15; /* 1/1000e */

 38: #define RE_CUT 3.
 39: /* < v, u_re * v * q > */
 40: static void f0_j_re(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 41:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 42:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 43:                     PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
 44: {
 45:   PetscReal n_e = PetscRealPart(u[0]);
 46:   if (dim==2) {
 47:     if (x[1] > RE_CUT || x[1] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
 48:       *f0 = n_e * 2.*PETSC_PI*x[0] * x[1] * constants[0]; /* n * r * v_|| * q */
 49:     } else {
 50:       *f0 = 0;
 51:     }
 52:   } else {
 53:     if (x[2] > RE_CUT || x[2] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
 54:       *f0 = n_e                * x[2] * constants[0];
 55:     } else {
 56:       *f0 = 0;
 57:     }
 58:   }
 59: }

 61: /* sum < v, u*v*q > */
 62: static void f0_jz_sum(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 63:                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 64:                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 65:                       PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar q[], PetscScalar *f0)
 66: {
 67:   PetscInt ii;
 68:   f0[0] = 0;
 69:   if (dim==2) {
 70:     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 */
 71:   } else {
 72:     for (ii=0;ii<Nf;ii++) f0[0] += u[ii]                * x[2] * q[ii]; /* n * v_|| * q  * v_0 */
 73:   }
 74: }

 76: /* < v, n_e > */
 77: static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 78:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 79:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 80:                   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 {
 85:     f0[0] =                        u[ii];
 86:   }
 87: }

 89: /* < v, n_e v_|| > */
 90: static void f0_vz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 91:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 92:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 93:                    PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
 94: {
 95:   PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
 96:   if (dim==2) f0[0] = u[ii] * 2.*PETSC_PI*x[0] * x[1]; /* n r v_|| */
 97:   else {
 98:     f0[0] =           u[ii] *                x[2]; /* n v_|| */
 99:   }
100: }

102: /* < v, n_e (v-shift) > */
103: static void f0_ve_shift(PetscInt dim, PetscInt Nf, PetscInt NfAux,
104:                          const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
105:                          const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
106:                          PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
107: {
108:   PetscReal vz = numConstants>0 ? PetscRealPart(constants[0]) : 0;
109:   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 */
110:   else {
111:     *f0 =           u[0] *                PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + (x[2]-vz)*(x[2]-vz)); /* n v */
112:   }
113: }

115:  /* CalculateE - Calculate the electric field  */
116:  /*  T        -- Electron temperature  */
117:  /*  n        -- Electron density  */
118:  /*  lnLambda --   */
119:  /*  eps0     --  */
120:  /*  E        -- output E, input \hat E */
121: static PetscReal CalculateE(PetscReal Tev, PetscReal n, PetscReal lnLambda, PetscReal eps0, PetscReal *E)
122: {
123:   PetscReal c,e,m;

126:   c = 299792458.0;
127:   e = 1.602176e-19;
128:   m = 9.10938e-31;
129:   if (1) {
130:     double Ec, Ehat = *E, betath = PetscSqrtReal(2*Tev*e/(m*c*c)), j0 = Ehat * 7/(PetscSqrtReal(2)*2) * PetscPowReal(betath,3) * n * e * c;
131:     Ec = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*c*c);
132:     *E = Ec;
133:     PetscPrintf(PETSC_COMM_WORLD, "CalculateE j0=%g Ec = %g\n",j0,Ec);
134:   } else {
135:     PetscReal Ed, vth;
136:     vth = PetscSqrtReal(8*Tev*e/(m*PETSC_PI));
137:     Ed =  n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*vth*vth);
138:     *E = Ed;
139:   }
140:   return(0);
141: }

143: static PetscReal Spitzer(PetscReal m_e, PetscReal e, PetscReal Z, PetscReal epsilon0,  PetscReal lnLam, PetscReal kTe_joules)
144: {
145:   PetscReal Fz = (1+1.198*Z+0.222*Z*Z)/(1+2.966*Z+0.753*Z*Z), eta;
146:   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);
147:   return eta;
148: }

150: /*  */
151: static PetscErrorCode testNone(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
152: {
154:   return(0);
155: }

157: /*  */
158: static PetscErrorCode testSpitzer(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
159: {
160:   PetscErrorCode    ierr;
161:   PetscInt          ii,nDMs;
162:   PetscDS           prob;
163:   static PetscReal  old_ratio = 1e10;
164:   TSConvergedReason reason;
165:   PetscReal         J,J_re,spit_eta,Te_kev=0,E,ratio,Z,n_e,v,v2;
166:   PetscScalar       user[2] = {0.,ctx->charges[0]}, q[LANDAU_MAX_SPECIES],tt[LANDAU_MAX_SPECIES],vz;
167:   PetscReal         dt;
168:   DM                pack, plexe = ctx->plex[0], plexi = (ctx->num_grids==1) ? NULL : ctx->plex[1];
169:   Vec               *XsubArray;

172:   if (ctx->num_species!=2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "ctx->num_species %D != 2",ctx->num_species);
173:   VecGetDM(X, &pack);
174:   if (!pack) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no DM");
175:   DMCompositeGetNumberDM(pack,&nDMs);
176:   if (nDMs != ctx->num_grids*ctx->batch_sz)  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nDMs != ctx->num_grids*ctx->batch_sz %D != %D",nDMs,ctx->num_grids*ctx->batch_sz);
177:   PetscMalloc(sizeof(*XsubArray)*nDMs, &XsubArray);
178:   DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
179:   TSGetTimeStep(ts,&dt);
180:   /* get current for each grid */
181:   for (ii=0;ii<ctx->num_species;ii++) q[ii] = ctx->charges[ii];
182:   DMGetDS(plexe, &prob);
183:   PetscDSSetConstants(prob, 2, &q[0]);
184:   PetscDSSetObjective(prob, 0, &f0_jz_sum);
185:   DMPlexComputeIntegralFEM(plexe,XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx,0) ],tt,NULL);
186:   J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
187:   if (plexi) { // add first (only) ion
188:     DMGetDS(plexi, &prob);
189:     PetscDSSetConstants(prob, 1, &q[1]);
190:     PetscDSSetObjective(prob, 0, &f0_jz_sum);
191:     DMPlexComputeIntegralFEM(plexi,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,1)],tt,NULL);
192:     J += -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
193:   }
194:   /* get N_e */
195:   DMGetDS(plexe, &prob);
196:   PetscDSSetConstants(prob, 1, user);
197:   PetscDSSetObjective(prob, 0, &f0_n);
198:   DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
199:   n_e = PetscRealPart(tt[0])*ctx->n_0;
200:   /* Z */
201:   Z = -ctx->charges[1]/ctx->charges[0];
202:   /* remove drift */
203:   if (0) {
204:     user[0] = 0; // electrons
205:     DMGetDS(plexe, &prob);
206:     PetscDSSetConstants(prob, 1, user);
207:     PetscDSSetObjective(prob, 0, &f0_vz);
208:     DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
209:     vz = ctx->n_0*PetscRealPart(tt[0])/n_e; /* non-dimensional */
210:   } else vz = 0;
211:   /* thermal velocity */
212:   DMGetDS(plexe, &prob);
213:   PetscDSSetConstants(prob, 1, &vz);
214:   PetscDSSetObjective(prob, 0, &f0_ve_shift);
215:   DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
216:   v = ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/n_e;   /* remove number density to get velocity */
217:   v2 = PetscSqr(v);                                    /* use real space: m^2 / s^2 */
218:   Te_kev = (v2*ctx->masses[0]*PETSC_PI/8)*kev_joul;    /* temperature in kev */
219:   //Te_kev = ctx->thermal_temps[0]/1.1604525e7;
220:   spit_eta = Spitzer(ctx->masses[0],-ctx->charges[0],Z,ctx->epsilon0,ctx->lnLam,Te_kev/kev_joul); /* kev --> J (kT) */
221:   if (0) {
222:     DMGetDS(plexe, &prob);
223:     PetscDSSetConstants(prob, 1, q);
224:     PetscDSSetObjective(prob, 0, &f0_j_re);
225:     DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
226:   } else tt[0] = 0;
227:   J_re = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
228:   DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
229:   PetscFree(XsubArray);

231:   if (rectx->use_spitzer_eta) {
232:     E = ctx->Ez = spit_eta*(rectx->j-J_re);
233:   } else {
234:     E = ctx->Ez; /* keep real E */
235:     rectx->j = J; /* cache */
236:   }

238:   ratio = E/J/spit_eta;
239:   if (stepi>10 && !rectx->use_spitzer_eta && (
240:         //(old_ratio-ratio < 1.e-3 && ratio > 0.99 && ratio < 1.01) ||
241:         //(old_ratio-ratio < 1.e-4 && ratio > 0.98 && ratio < 1.02) ||
242:         (old_ratio-ratio < 1.e-6))) {
243:     rectx->pulse_start = time + 0.98*dt;
244:     rectx->use_spitzer_eta = PETSC_TRUE;
245:   }
246:   TSGetConvergedReason(ts,&reason);
247:   TSGetConvergedReason(ts,&reason);
248:   if ((rectx->plotting) || stepi == 0 || reason || rectx->pulse_start == time + 0.98*dt) {
249:     PetscPrintf(ctx->comm, "testSpitzer: %4D) 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,time,n_e/ctx->n_0,ctx->Ez,J,J_re,100*J_re/J, Te_kev,Z,ratio,old_ratio-ratio, rectx->use_spitzer_eta ? "using Spitzer eta*J E" : "constant E",rectx->pulse_start != time + 0.98*dt ? "normal" : "transition",spit_eta);
250:     if (rectx->pulse_start == time + 0.98*dt) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spitzer complete ratio=%g",ratio);
251:   }
252:   old_ratio = ratio;
253:   return(0);
254: }

256: static const double ppp = 2;
257: static void f0_0_diff_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
258:                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
259:                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
260:                           PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
261: {
262:   LandauCtx       *ctx = (LandauCtx *)constants;
263:   REctx           *rectx = (REctx*)ctx->data;
264:   PetscInt        ii = rectx->idx, i;
265:   const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */
266:   const PetscReal n = ctx->n[ii];
267:   PetscReal       diff, f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
268:   for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
269:   f_maxwell = n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
270:   diff = 2.*PETSC_PI*x[0]*(PetscRealPart(u[ii]) - f_maxwell);
271:   f0[0] = PetscPowReal(diff,ppp);
272: }
273: static void f0_0_maxwellian_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
274:                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
275:                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
276:                           PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
277: {
278:   LandauCtx       *ctx = (LandauCtx *)constants;
279:   REctx           *rectx = (REctx*)ctx->data;
280:   PetscInt        ii = rectx->idx, i;
281:   const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */
282:   const PetscReal n = ctx->n[ii];
283:   PetscReal       f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
284:   for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
285:   f_maxwell = 2.*PETSC_PI*x[0] * n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
286:   f0[0] = PetscPowReal(f_maxwell,ppp);
287: }

289: /*  */
290: static PetscErrorCode testStable(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
291: {
292:   PetscErrorCode    ierr;
293:   PetscDS           prob;
294:   Vec               X2;
295:   PetscReal         ediff,idiff=0,lpm0,lpm1=1;
296:   PetscScalar       tt[LANDAU_MAX_SPECIES];
297:   DM                dm, plex = ctx->plex[0];

300:   VecGetDM(X, &dm);
301:   DMGetDS(plex, &prob);
302:   VecDuplicate(X,&X2);
303:   VecCopy(X,X2);
304:   if (!rectx->X_0) {
305:     VecDuplicate(X,&rectx->X_0);
306:     VecCopy(X,rectx->X_0);
307:   }
308:   VecAXPY(X,-1.0,rectx->X_0);
309:   PetscDSSetConstants(prob, sizeof(LandauCtx)/sizeof(PetscScalar), (PetscScalar*)ctx);
310:   rectx->idx = 0;
311:   PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
312:   DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
313:   ediff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
314:   PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
315:   DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
316:   lpm0 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
317:   if (ctx->num_species>1) {
318:     rectx->idx = 1;
319:     PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
320:     DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
321:     idiff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
322:     PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
323:     DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
324:     lpm1 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
325:   }
326:   PetscPrintf(PETSC_COMM_WORLD, "%s %D) time=%10.3e n-%d norm electrons/max=%20.13e ions/max=%20.13e\n", "----",stepi,time,(int)ppp,ediff/lpm0,idiff/lpm1);
327:   /* view */
328:   VecCopy(X2,X);
329:   VecDestroy(&X2);
330:   if (islast) {
331:     VecDestroy(&rectx->X_0);
332:     rectx->X_0 = NULL;
333:   }
334:   return(0);
335: }

337: static PetscErrorCode EInduction(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
338: {
339:   REctx             *rectx = (REctx*)ctx->data;
340:   PetscErrorCode    ierr;
341:   PetscInt          ii;
342:   DM                dm,plex;
343:   PetscScalar       tt[LANDAU_MAX_SPECIES], qv0[LANDAU_MAX_SPECIES];
344:   PetscReal         dJ_dt;
345:   PetscDS           prob;

348:   for (ii=0;ii<ctx->num_species;ii++) qv0[ii] = ctx->charges[ii]*ctx->v_0;
349:   VecGetDM(X, &dm);
350:   DMGetDS(dm, &prob);
351:   DMConvert(dm, DMPLEX, &plex);
352:   /* get d current / dt */
353:   PetscDSSetConstants(prob, ctx->num_species, qv0);
354:   PetscDSSetObjective(prob, 0, &f0_jz_sum);
355:   if (!X_t) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "X_t");
356:   DMPlexComputeIntegralFEM(plex,X_t,tt,NULL);
357:   dJ_dt = -ctx->n_0*PetscRealPart(tt[0])/ctx->t_0;
358:   /* E induction */
359:   *a_E = -rectx->L*dJ_dt + rectx->Ez_initial;
360:   DMDestroy(&plex);
361:   return(0);
362: }

364: static PetscErrorCode EConstant(Vec X,  Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
365: {
367:   *a_E = ctx->Ez;
368:   return(0);
369: }

371: static PetscErrorCode ENone(Vec X,  Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
372: {
374:   *a_E = 0;
375:   return(0);
376: }

378: /* ------------------------------------------------------------------- */
379: /*
380:    FormSource - Evaluates source terms F(t).

382:    Input Parameters:
383: .  ts - the TS context
384: .  time -
385: .  X_dummmy - input vector
386: .  dummy - optional user-defined context, as set by SNESSetFunction()

388:    Output Parameter:
389: .  F - function vector
390:  */
391: static PetscErrorCode FormSource(TS ts, PetscReal ftime, Vec X_dummmy, Vec F, void *dummy)
392: {
393:   PetscReal      new_imp_rate;
394:   LandauCtx      *ctx;
395:   DM             pack;
397:   REctx          *rectx;

400:   TSGetDM(ts,&pack);
401:   DMGetApplicationContext(pack, &ctx);
402:   rectx = (REctx*)ctx->data;
403:   /* check for impurities */
404:   rectx->impuritySrcRate(ftime,&new_imp_rate,ctx);
405:   if (new_imp_rate != 0) {
406:     if (new_imp_rate != rectx->current_rate) {
407:       PetscInt       ii;
408:       PetscReal      dne_dt,dni_dt,tilda_ns[LANDAU_MAX_SPECIES],temps[LANDAU_MAX_SPECIES];
409:       Vec            globFarray[LANDAU_MAX_GRIDS*LANDAU_MAX_BATCH_SZ];
410:       rectx->current_rate = new_imp_rate;
411:       for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) tilda_ns[ii] = 0;
412:       for (ii=1;ii<LANDAU_MAX_SPECIES;ii++)    temps[ii] = 1;
413:       dni_dt = new_imp_rate               /* *ctx->t_0 */; /* fully ionized immediately, no normalize, stay in non-dim */
414:       dne_dt = new_imp_rate*rectx->Ne_ion /* *ctx->t_0 */;
415:       tilda_ns[0] = dne_dt;        tilda_ns[rectx->imp_idx] = dni_dt;
416:       temps[0]    = rectx->T_cold;    temps[rectx->imp_idx] = rectx->T_cold;
417:       PetscInfo4(ctx->plex[0], "\tHave new_imp_rate= %10.3e time= %10.3e de/dt= %10.3e di/dt= %10.3e ***\n",new_imp_rate,ftime,dne_dt,dni_dt);
418:       DMCompositeGetAccessArray(pack, F, ctx->num_grids*ctx->batch_sz, NULL, globFarray);
419:       for (PetscInt grid=0 ; grid<ctx->num_grids ; grid++) {
420:         /* add it */
421:         LandauAddMaxwellians(ctx->plex[grid],globFarray[ LAND_PACK_IDX(0,grid) ],ftime,temps,tilda_ns,grid,0,ctx);
422:         VecViewFromOptions(globFarray[ LAND_PACK_IDX(0,grid) ],NULL,"-vec_view_sources");
423:       }
424:       // Does DMCompositeRestoreAccessArray copy the data back? (no)
425:       DMCompositeRestoreAccessArray(pack, F, ctx->num_grids*ctx->batch_sz, NULL, globFarray);
426:     }
427:   } else {
428:     VecZeroEntries(F);
429:     rectx->current_rate = 0;
430:   }
431:   return(0);
432: }
433: PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx)
434: {
435:   LandauCtx         *ctx = (LandauCtx*) actx;   /* user-defined application context */
436:   REctx             *rectx = (REctx*)ctx->data;
437:   DM                pack;
438:   Vec               globXArray[LANDAU_MAX_GRIDS*LANDAU_MAX_BATCH_SZ];
439:   TSConvergedReason reason;
440:   PetscErrorCode    ierr;
442:   VecGetDM(X, &pack);
443:   DMCompositeGetAccessArray(pack, X, ctx->num_grids*ctx->batch_sz, NULL, globXArray);
444:   if (stepi > rectx->plotStep && rectx->plotting) {
445:     rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */
446:     rectx->plotIdx++;
447:   }
448:   /* view */
449:   TSGetConvergedReason(ts,&reason);
450:   if (time/rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) {
451:     if ((reason || stepi==0 || rectx->plotIdx%rectx->print_period==0) && ctx->verbose > 0) {
452:       /* print norms */
453:       LandauPrintNorms(X, stepi);
454:     }
455:     if (!rectx->plotting) { /* first step of possible backtracks */
456:       rectx->plotting = PETSC_TRUE;
457:       /* diagnostics + change E field with Sptizer (not just a monitor) */
458:       rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
459:     } else {
460:       PetscPrintf(PETSC_COMM_WORLD, "\t\t ERROR SKIP test spit ------\n");
461:       rectx->plotting = PETSC_TRUE;
462:     }
463:     PetscObjectSetName((PetscObject) globXArray[ LAND_PACK_IDX(ctx->batch_view_idx,rectx->grid_view_idx) ], rectx->grid_view_idx==0 ? "ue" : "ui");
464:     /* view, overwrite step when back tracked */
465:     DMSetOutputSequenceNumber(pack, rectx->plotIdx, time*ctx->t_0);
466:     VecViewFromOptions(globXArray[ LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx) ],NULL,"-vec_view");

468:     rectx->plotStep = stepi;
469:   } else {
470:     if (rectx->plotting) PetscPrintf(PETSC_COMM_WORLD," ERROR rectx->plotting=%D step %D\n",rectx->plotting,stepi);
471:     /* diagnostics + change E field with Sptizer (not just a monitor) - can we lag this? */
472:     rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
473:   }
474:   /* parallel check that only works of all batches are identical */
475:   if (reason && ctx->verbose > 3) {
476:     PetscReal    val,rval;
477:     PetscMPIInt  rank;
478:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
479:     for (PetscInt grid=0;grid<ctx->num_grids;grid++) {
480:       PetscInt nerrors=0;
481:       for (PetscInt i=0; i<ctx->batch_sz;i++) {
482:         VecNorm(globXArray[ LAND_PACK_IDX(i,grid) ],NORM_2,&val);
483:         if (i==0) rval = val;
484:         else if ((val=PetscAbs(val-rval)/rval) > 1000*PETSC_MACHINE_EPSILON) {
485:           PetscPrintf(PETSC_COMM_SELF, " [%D] Warning %D.%D) diff = %2.15e\n",rank,grid,i,val);
486:           nerrors++;
487:         }
488:       }
489:       if (nerrors) {
490:         PetscPrintf(PETSC_COMM_SELF, " ***** [%D] ERROR max %D errors\n",rank,nerrors);
491:       } else {
492:         PetscPrintf(PETSC_COMM_WORLD, "[%D] %D) batch consistency check OK\n",rank,grid);
493:       }
494:     }
495:   }
496:   rectx->idx = 0;
497:   DMCompositeRestoreAccessArray(pack, X, ctx->num_grids*ctx->batch_sz, NULL, globXArray);
498:   return(0);
499: }

501: PetscErrorCode PreStep(TS ts)
502: {
504:   LandauCtx      *ctx;
505:   REctx          *rectx;
506:   DM             dm;
507:   PetscInt       stepi;
508:   PetscReal      time;
509:   Vec            X;

512:   /* not used */
513:   TSGetDM(ts,&dm);
514:   TSGetTime(ts,&time);
515:   TSGetSolution(ts,&X);
516:   DMGetApplicationContext(dm, &ctx);
517:   rectx = (REctx*)ctx->data;
518:   TSGetStepNumber(ts, &stepi);
519:   /* update E */
520:   rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez);
521:   return(0);
522: }

524: /* 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) */
525: static PetscErrorCode stepSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
526: {
527:   REctx         *rectx = (REctx*)ctx->data;

530:   if (time >= rectx->pulse_start) *rho = rectx->pulse_rate;
531:   else *rho = 0.;
532:   return(0);
533: }
534: static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
535: {
537:   *rho = 0.;
538:   return(0);
539: }
540: static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
541: {
542:   REctx *rectx = (REctx*)ctx->data;

545:   if (rectx->pulse_start == PETSC_MAX_REAL) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"'-ex2_pulse_start_time X' must be used with '-ex2_impurity_source_type pulse'");
546:   if (time < rectx->pulse_start || time > rectx->pulse_start + 3*rectx->pulse_width) *rho = 0;
547:   /* else if (0) { */
548:   /*   double t = time - rectx->pulse_start, start = rectx->pulse_width, stop = 2*rectx->pulse_width, cycle = 3*rectx->pulse_width, steep = 5, xi = 0.75 - (stop - start)/(2* cycle); */
549:   /*   *rho = rectx->pulse_rate * (cycle / (stop - start)) / (1 + PetscExpReal(steep*(PetscSinReal(2*PETSC_PI*((t - start)/cycle + xi)) - PetscSinReal(2*PETSC_PI*xi)))); */
550:   /* } else if (0) { */
551:   /*   double x = 2*(time - rectx->pulse_start)/(3*rectx->pulse_width) - 1; */
552:   /*   if (x==1 || x==-1) *rho = 0; */
553:   /*   else *rho = rectx->pulse_rate * PetscExpReal(-1/(1-x*x)); */
554:   /* } */
555:   else {
556:     double x = PetscSinReal((time-rectx->pulse_start)/(3*rectx->pulse_width)*2*PETSC_PI - PETSC_PI/2) + 1; /* 0:2, integrates to 1.0 */
557:     *rho = rectx->pulse_rate * x / (3*rectx->pulse_width);
558:     if (!rectx->use_spitzer_eta) rectx->use_spitzer_eta = PETSC_TRUE; /* use it next time */
559:   }
560:   return(0);
561: }

565: static PetscErrorCode ProcessREOptions(REctx *rectx, const LandauCtx *ctx, DM dm, const char prefix[])
566: {
567:   PetscErrorCode    ierr;
568:   PetscFunctionList plist = NULL, testlist = NULL, elist = NULL;
569:   char              pname[256],testname[256],ename[256];
570:   DM                dm_dummy;
571:   PetscBool         Connor_E = PETSC_FALSE;

574:   DMCreate(PETSC_COMM_WORLD,&dm_dummy);
575:   rectx->Ne_ion = 1;                 /* number of electrons given up by impurity ion */
576:   rectx->T_cold = .005;              /* kev */
577:   rectx->ion_potential = 15;         /* ev */
578:   rectx->L = 2;
579:   rectx->X_0 = NULL;
580:   rectx->imp_idx = ctx->num_species - 1; /* default ionized impurity as last one */
581:   rectx->pulse_start = PETSC_MAX_REAL;
582:   rectx->pulse_width = 1;
583:   rectx->plotStep = PETSC_MAX_INT;
584:   rectx->pulse_rate = 1.e-1;
585:   rectx->current_rate = 0;
586:   rectx->plotIdx = 0;
587:   rectx->j = 0;
588:   rectx->plotDt = 1.0;
589:   rectx->plotting = PETSC_FALSE;
590:   rectx->use_spitzer_eta = PETSC_FALSE;
591:   rectx->idx = 0;
592:   rectx->print_period = 10;
593:   rectx->grid_view_idx = 0;
594:   /* Register the available impurity sources */
595:   PetscFunctionListAdd(&plist,"step",&stepSrc);
596:   PetscFunctionListAdd(&plist,"none",&zeroSrc);
597:   PetscFunctionListAdd(&plist,"pulse",&pulseSrc);
598:   PetscStrcpy(pname,"none");
599:   PetscFunctionListAdd(&testlist,"none",&testNone);
600:   PetscFunctionListAdd(&testlist,"spitzer",&testSpitzer);
601:   PetscFunctionListAdd(&testlist,"stable",&testStable);
602:   PetscStrcpy(testname,"none");
603:   PetscFunctionListAdd(&elist,"none",&ENone);
604:   PetscFunctionListAdd(&elist,"induction",&EInduction);
605:   PetscFunctionListAdd(&elist,"constant",&EConstant);
606:   PetscStrcpy(ename,"constant");

608:   PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");
609:   PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "ex2.c", rectx->plotDt, &rectx->plotDt, NULL);
610:   if (rectx->plotDt < 0) rectx->plotDt = 1e30;
611:   if (rectx->plotDt == 0) rectx->plotDt = 1e-30;
612:   PetscOptionsInt("-ex2_print_period", "Plotting interval", "ex2.c", rectx->print_period, &rectx->print_period, NULL);
613:   PetscOptionsInt("-ex2_grid_view_idx", "grid_view_idx", "ex2.c", rectx->grid_view_idx, &rectx->grid_view_idx, NULL);
614:   if (rectx->grid_view_idx >= ctx->num_grids) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"rectx->grid_view_idx (%D) >= ctx->num_grids (%D)",rectx->imp_idx,ctx->num_grids);
615:   PetscOptionsFList("-ex2_impurity_source_type","Name of impurity source to run","",plist,pname,pname,sizeof(pname),NULL);
616:   PetscOptionsFList("-ex2_test_type","Name of test to run","",testlist,testname,testname,sizeof(testname),NULL);
617:   PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL);
618:   if ((rectx->imp_idx >= ctx->num_species || rectx->imp_idx < 1) && ctx->num_species > 1) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"index of sink for impurities ions is out of range (%D), must be > 0 && < NS",rectx->imp_idx);
619:   PetscOptionsFList("-ex2_e_field_type","Electric field type","",elist,ename,ename,sizeof(ename),NULL);
620:   rectx->Ne_ion = -ctx->charges[rectx->imp_idx]/ctx->charges[0];
621:   PetscOptionsReal("-ex2_t_cold","Temperature of cold electron and ions after ionization in keV","none",rectx->T_cold,&rectx->T_cold, NULL);
622:   PetscOptionsReal("-ex2_pulse_start_time","Time at which pulse happens for 'pulse' source","none",rectx->pulse_start,&rectx->pulse_start, NULL);
623:   PetscOptionsReal("-ex2_pulse_width_time","Width of pulse 'pulse' source","none",rectx->pulse_width,&rectx->pulse_width, NULL);
624:   PetscOptionsReal("-ex2_pulse_rate","Number density of pulse for 'pulse' source","none",rectx->pulse_rate,&rectx->pulse_rate, NULL);
625:   rectx->T_cold *= 1.16e7; /* convert to Kelvin */
626:   PetscOptionsReal("-ex2_ion_potential","Potential to ionize impurity (should be array) in ev","none",rectx->ion_potential,&rectx->ion_potential, NULL);
627:   PetscOptionsReal("-ex2_inductance","Inductance E feild","none",rectx->L,&rectx->L, NULL);
628:   PetscOptionsBool("-ex2_connor_e_field_units","Scale Ex but Connor-Hastie E_c","none",Connor_E,&Connor_E, NULL);
629:   PetscInfo5(dm_dummy, "Num electrons from ions=%g, T_cold=%10.3e, ion potential=%10.3e, E_z=%10.3e v_0=%10.3e\n",rectx->Ne_ion,rectx->T_cold,rectx->ion_potential,ctx->Ez,ctx->v_0);
630:   PetscOptionsEnd();
631:   /* get impurity source rate function */
632:   PetscFunctionListFind(plist,pname,&rectx->impuritySrcRate);
633:   if (!rectx->impuritySrcRate) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No impurity source function found '%s'",pname);
634:   PetscFunctionListFind(testlist,testname,&rectx->test);
635:   if (!rectx->test) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No test found '%s'",testname);
636:   PetscFunctionListFind(elist,ename,&rectx->E);
637:   if (!rectx->E) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No E field function found '%s'",ename);
638:   PetscFunctionListDestroy(&plist);
639:   PetscFunctionListDestroy(&testlist);
640:   PetscFunctionListDestroy(&elist);

642:   /* convert E from Connor-Hastie E_c units to real if doing Spitzer E */
643:   if (Connor_E) {
644:     PetscReal E = ctx->Ez, Tev = ctx->thermal_temps[0]*8.621738e-5, n = ctx->n_0*ctx->n[0];
645:     CalculateE(Tev, n, ctx->lnLam, ctx->epsilon0, &E);
646:     ((LandauCtx*)ctx)->Ez *= E;
647:   }
648:   DMDestroy(&dm_dummy);
649:   return(0);
650: }

652: static PetscErrorCode MatrixNfDestroy(void *ptr)
653: {
654:   PetscInt *nf = (PetscInt *)ptr;
655:   PetscErrorCode  ierr;
657:   PetscFree(nf);
658:   return(0);
659: }

661: int main(int argc, char **argv)
662: {
663:   DM             pack;
664:   Vec            X,*XsubArray;
666:   PetscInt       dim = 2, nDMs;
667:   TS             ts;
668:   Mat            J;
669:   PetscDS        prob;
670:   LandauCtx      *ctx;
671:   REctx          *rectx;
672: #if defined PETSC_USE_LOG
673:   PetscLogStage  stage;
674: #endif
675:   PetscMPIInt    rank;
676: #if defined(PETSC_HAVE_THREADSAFETY)
677:   double         starttime, endtime;
678: #endif
679:   PetscContainer container;
680:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
681:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
682:   if (rank) { /* turn off output stuff for duplicate runs */
683:     PetscOptionsClearValue(NULL,"-dm_view");
684:     PetscOptionsClearValue(NULL,"-vec_view");
685:     PetscOptionsClearValue(NULL,"-dm_view_diff");
686:     PetscOptionsClearValue(NULL,"-vec_view_diff");
687:     PetscOptionsClearValue(NULL,"-dm_view_sources");
688:     PetscOptionsClearValue(NULL,"-vec_view_0");
689:     PetscOptionsClearValue(NULL,"-dm_view_0");
690:     PetscOptionsClearValue(NULL,"-vec_view_sources");
691:     PetscOptionsClearValue(NULL,"-info"); /* this does not work */
692:   }
693:   PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);
694:   /* Create a mesh */
695:   LandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, "", &X, &J, &pack);
696:   DMCompositeGetNumberDM(pack,&nDMs);
697:   PetscMalloc(sizeof(*XsubArray)*nDMs, &XsubArray);
698:   PetscObjectSetName((PetscObject)J, "Jacobian");
699:   PetscObjectSetName((PetscObject)X, "f");
700:   LandauCreateMassMatrix(pack, NULL);
701:   DMGetApplicationContext(pack, &ctx);
702:   DMSetUp(pack);
703:   /* context */
704:   PetscNew(&rectx);
705:   ctx->data = rectx;
706:   ProcessREOptions(rectx,ctx,pack,"");
707:   DMGetDS(pack, &prob);
708:   DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
709:   PetscObjectSetName((PetscObject) XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx) ], rectx->grid_view_idx==0 ? "ue" : "ui");
710:   DMViewFromOptions(ctx->plex[rectx->grid_view_idx],NULL,"-dm_view");
711:   DMViewFromOptions(ctx->plex[rectx->grid_view_idx], NULL,"-dm_view_0");
712:   VecViewFromOptions(XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx,rectx->grid_view_idx) ], NULL,"-vec_view_0"); // initial condition (monitor plots after step)
713:   DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
714:   PetscFree(XsubArray);
715:   VecViewFromOptions(X, NULL,"-vec_view_global"); // initial condition (monitor plots after step)
716:   DMSetOutputSequenceNumber(pack, 0, 0.0);
717:   /* Create timestepping solver context */
718:   TSCreate(PETSC_COMM_SELF,&ts);
719:   TSSetDM(ts,pack);
720:   TSSetIFunction(ts,NULL,LandauIFunction,NULL);
721:   TSSetIJacobian(ts,J,J,LandauIJacobian,NULL);
722:   TSSetRHSFunction(ts,NULL,FormSource,NULL);
723:   TSSetFromOptions(ts);
724:   TSSetSolution(ts,X);
725:   TSSetApplicationContext(ts, ctx);
726:   TSMonitorSet(ts,Monitor,ctx,NULL);
727:   TSSetPreStep(ts,PreStep);
728:   {
729:     PetscInt *pNf;
730:     PetscContainerCreate(PETSC_COMM_SELF, &container);
731:     PetscMalloc(sizeof(*pNf), &pNf);
732:     *pNf = ctx->batch_sz + 100000*ctx->numConcurrency;
733:     PetscContainerSetPointer(container, (void *)pNf);
734:     PetscContainerSetUserDestroy(container, MatrixNfDestroy);
735:     PetscObjectCompose((PetscObject)ctx->J, "batch size", (PetscObject) container);
736:     PetscContainerDestroy(&container);
737:   }
738:   rectx->Ez_initial = ctx->Ez;       /* cache for induction caclulation - applied E field */
739:   if (1) { /* warm up an test just LandauIJacobian */
740:     Vec           vec;
741:     PetscInt      nsteps;
742:     PetscReal     dt;
743:     PetscLogStageRegister("Warmup", &stage);
744:     PetscLogStagePush(stage);
745:     VecDuplicate(X,&vec);
746:     VecCopy(X,vec);
747:     TSGetMaxSteps(ts,&nsteps);
748:     TSGetTimeStep(ts,&dt);
749:     TSSetMaxSteps(ts,1);
750:     TSSolve(ts,X);
751:     TSSetMaxSteps(ts,nsteps);
752:     TSSetStepNumber(ts,0);
753:     TSSetTime(ts,0);
754:     TSSetTimeStep(ts,dt);
755:     rectx->plotIdx = 0;
756:     rectx->plotting = PETSC_FALSE;
757:     PetscLogStagePop();
758:     VecCopy(vec,X);
759:     VecDestroy(&vec);
760:     ctx->aux_bool = PETSC_FALSE; // flag for not a clean Jacobian
761:   }
762:   /* go */
763:   PetscLogStageRegister("Solve", &stage);
764:   ctx->stage = 0; // lets not use this stage
765: #if defined(PETSC_HAVE_THREADSAFETY)
766:   ctx->stage = 1; // not set with thread safty
767: #endif
768:   TSSetSolution(ts,X);
769:   PetscLogStagePush(stage);
770: #if defined(PETSC_HAVE_THREADSAFETY)
771:   starttime = MPI_Wtime();
772: #endif
773:   TSSolve(ts,X);
774:   PetscLogStagePop();
775: #if defined(PETSC_HAVE_THREADSAFETY)
776:   endtime = MPI_Wtime();
777:   ctx->times[LANDAU_EX2_TSSOLVE] += (endtime - starttime);
778: #endif
779:   VecViewFromOptions(X, NULL,"-vec_view_global");
780:   /* clean up */
781:   LandauDestroyVelocitySpace(&pack);
782:   TSDestroy(&ts);
783:   VecDestroy(&X);
784:   PetscFree(rectx);
785:   PetscContainerDestroy(&container);
786:   PetscFinalize();
787:   return ierr;
788: }

790: /*TEST

792:   testset:
793:     requires: p4est !complex double
794:     output_file: output/ex2_0.out
795:     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-10 -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 -1 -ts_rtol 1e-3 -ts_dt 1.e-1 -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 -pc_type lu -ksp_type preonly -dm_landau_amr_levels_max 2,2 -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 1 -dm_landau_gpu_assembly true -dm_landau_batch_size 2
796:     test:
797:       suffix: cpu
798:       args: -dm_landau_device_type cpu
799:     test:
800:       suffix: kokkos
801:       requires: kokkos_kernels
802:       args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos
803:     test:
804:       suffix: cuda
805:       requires: cuda
806:       args: -dm_landau_device_type cuda -dm_mat_type aijcusparse -dm_vec_type cuda -mat_cusparse_use_cpu_solve

808: TEST*/