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;
125: c = 299792458.0;
126: e = 1.602176e-19;
127: m = 9.10938e-31;
128: if (1) {
129: double Ec, Ehat = *E, betath = PetscSqrtReal(2*Tev*e/(m*c*c)), j0 = Ehat * 7/(PetscSqrtReal(2)*2) * PetscPowReal(betath,3) * n * e * c;
130: Ec = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*c*c);
131: *E = Ec;
132: PetscPrintf(PETSC_COMM_WORLD, "CalculateE j0=%g Ec = %g\n",j0,Ec);
133: } else {
134: PetscReal Ed, vth;
135: vth = PetscSqrtReal(8*Tev*e/(m*PETSC_PI));
136: Ed = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*vth*vth);
137: *E = Ed;
138: }
139: return 0;
140: }
142: static PetscReal Spitzer(PetscReal m_e, PetscReal e, PetscReal Z, PetscReal epsilon0, PetscReal lnLam, PetscReal kTe_joules)
143: {
144: PetscReal Fz = (1+1.198*Z+0.222*Z*Z)/(1+2.966*Z+0.753*Z*Z), eta;
145: 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);
146: return eta;
147: }
149: /* */
150: static PetscErrorCode testNone(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
151: {
153: return 0;
154: }
156: /* */
157: static PetscErrorCode testSpitzer(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
158: {
159: PetscInt ii,nDMs;
160: PetscDS prob;
161: static PetscReal old_ratio = 1e10;
162: TSConvergedReason reason;
163: PetscReal J,J_re,spit_eta,Te_kev=0,E,ratio,Z,n_e,v,v2;
164: PetscScalar user[2] = {0.,ctx->charges[0]}, q[LANDAU_MAX_SPECIES],tt[LANDAU_MAX_SPECIES],vz;
165: PetscReal dt;
166: DM pack, plexe = ctx->plex[0], plexi = (ctx->num_grids==1) ? NULL : ctx->plex[1];
167: Vec *XsubArray;
171: VecGetDM(X, &pack);
173: DMCompositeGetNumberDM(pack,&nDMs);
175: PetscMalloc(sizeof(*XsubArray)*nDMs, &XsubArray);
176: DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
177: TSGetTimeStep(ts,&dt);
178: /* get current for each grid */
179: for (ii=0;ii<ctx->num_species;ii++) q[ii] = ctx->charges[ii];
180: DMGetDS(plexe, &prob);
181: PetscDSSetConstants(prob, 2, &q[0]);
182: PetscDSSetObjective(prob, 0, &f0_jz_sum);
183: DMPlexComputeIntegralFEM(plexe,XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx,0) ],tt,NULL);
184: J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
185: if (plexi) { // add first (only) ion
186: DMGetDS(plexi, &prob);
187: PetscDSSetConstants(prob, 1, &q[1]);
188: PetscDSSetObjective(prob, 0, &f0_jz_sum);
189: DMPlexComputeIntegralFEM(plexi,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,1)],tt,NULL);
190: J += -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
191: }
192: /* get N_e */
193: DMGetDS(plexe, &prob);
194: PetscDSSetConstants(prob, 1, user);
195: PetscDSSetObjective(prob, 0, &f0_n);
196: DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
197: n_e = PetscRealPart(tt[0])*ctx->n_0;
198: /* Z */
199: Z = -ctx->charges[1]/ctx->charges[0];
200: /* remove drift */
201: if (0) {
202: user[0] = 0; // electrons
203: DMGetDS(plexe, &prob);
204: PetscDSSetConstants(prob, 1, user);
205: PetscDSSetObjective(prob, 0, &f0_vz);
206: DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
207: vz = ctx->n_0*PetscRealPart(tt[0])/n_e; /* non-dimensional */
208: } else vz = 0;
209: /* thermal velocity */
210: DMGetDS(plexe, &prob);
211: PetscDSSetConstants(prob, 1, &vz);
212: PetscDSSetObjective(prob, 0, &f0_ve_shift);
213: DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
214: v = ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/n_e; /* remove number density to get velocity */
215: v2 = PetscSqr(v); /* use real space: m^2 / s^2 */
216: Te_kev = (v2*ctx->masses[0]*PETSC_PI/8)*kev_joul; /* temperature in kev */
217: spit_eta = Spitzer(ctx->masses[0],-ctx->charges[0],Z,ctx->epsilon0,ctx->lnLam,Te_kev/kev_joul); /* kev --> J (kT) */
218: if (0) {
219: DMGetDS(plexe, &prob);
220: PetscDSSetConstants(prob, 1, q);
221: PetscDSSetObjective(prob, 0, &f0_j_re);
222: DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);
223: } else tt[0] = 0;
224: J_re = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
225: DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
226: PetscFree(XsubArray);
228: if (rectx->use_spitzer_eta) {
229: E = ctx->Ez = spit_eta*(rectx->j-J_re);
230: } else {
231: E = ctx->Ez; /* keep real E */
232: rectx->j = J; /* cache */
233: }
235: ratio = E/J/spit_eta;
236: if (stepi>10 && !rectx->use_spitzer_eta && (
237: (old_ratio-ratio < 1.e-6))) {
238: rectx->pulse_start = time + 0.98*dt;
239: rectx->use_spitzer_eta = PETSC_TRUE;
240: }
241: TSGetConvergedReason(ts,&reason);
242: TSGetConvergedReason(ts,&reason);
243: if ((rectx->plotting) || stepi == 0 || reason || rectx->pulse_start == time + 0.98*dt) {
244: 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);
246: }
247: old_ratio = ratio;
248: return 0;
249: }
251: static const double ppp = 2;
252: static void f0_0_diff_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
253: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
254: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
255: 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 diff, 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 = n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
265: diff = 2.*PETSC_PI*x[0]*(PetscRealPart(u[ii]) - f_maxwell);
266: f0[0] = PetscPowReal(diff,ppp);
267: }
268: static void f0_0_maxwellian_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
269: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
270: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
271: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
272: {
273: LandauCtx *ctx = (LandauCtx *)constants;
274: REctx *rectx = (REctx*)ctx->data;
275: PetscInt ii = rectx->idx, i;
276: const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */
277: const PetscReal n = ctx->n[ii];
278: PetscReal f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
279: for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
280: f_maxwell = 2.*PETSC_PI*x[0] * n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
281: f0[0] = PetscPowReal(f_maxwell,ppp);
282: }
284: /* */
285: static PetscErrorCode testStable(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
286: {
287: PetscDS prob;
288: Vec X2;
289: PetscReal ediff,idiff=0,lpm0,lpm1=1;
290: PetscScalar tt[LANDAU_MAX_SPECIES];
291: DM dm, plex = ctx->plex[0];
294: VecGetDM(X, &dm);
295: DMGetDS(plex, &prob);
296: VecDuplicate(X,&X2);
297: VecCopy(X,X2);
298: if (!rectx->X_0) {
299: VecDuplicate(X,&rectx->X_0);
300: VecCopy(X,rectx->X_0);
301: }
302: VecAXPY(X,-1.0,rectx->X_0);
303: PetscDSSetConstants(prob, sizeof(LandauCtx)/sizeof(PetscScalar), (PetscScalar*)ctx);
304: rectx->idx = 0;
305: PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
306: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
307: ediff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
308: PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
309: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
310: lpm0 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
311: if (ctx->num_species>1) {
312: rectx->idx = 1;
313: PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
314: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
315: idiff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
316: PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
317: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
318: lpm1 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
319: }
320: 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);
321: /* view */
322: VecCopy(X2,X);
323: VecDestroy(&X2);
324: if (islast) {
325: VecDestroy(&rectx->X_0);
326: rectx->X_0 = NULL;
327: }
328: return 0;
329: }
331: static PetscErrorCode EInduction(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
332: {
333: REctx *rectx = (REctx*)ctx->data;
334: PetscInt ii;
335: DM dm,plex;
336: PetscScalar tt[LANDAU_MAX_SPECIES], qv0[LANDAU_MAX_SPECIES];
337: PetscReal dJ_dt;
338: PetscDS prob;
341: for (ii=0;ii<ctx->num_species;ii++) qv0[ii] = ctx->charges[ii]*ctx->v_0;
342: VecGetDM(X, &dm);
343: DMGetDS(dm, &prob);
344: DMConvert(dm, DMPLEX, &plex);
345: /* get d current / dt */
346: PetscDSSetConstants(prob, ctx->num_species, qv0);
347: PetscDSSetObjective(prob, 0, &f0_jz_sum);
349: DMPlexComputeIntegralFEM(plex,X_t,tt,NULL);
350: dJ_dt = -ctx->n_0*PetscRealPart(tt[0])/ctx->t_0;
351: /* E induction */
352: *a_E = -rectx->L*dJ_dt + rectx->Ez_initial;
353: DMDestroy(&plex);
354: return 0;
355: }
357: static PetscErrorCode EConstant(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
358: {
360: *a_E = ctx->Ez;
361: return 0;
362: }
364: static PetscErrorCode ENone(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
365: {
367: *a_E = 0;
368: return 0;
369: }
371: /* ------------------------------------------------------------------- */
372: /*
373: FormSource - Evaluates source terms F(t).
375: Input Parameters:
376: . ts - the TS context
377: . time -
378: . X_dummmy - input vector
379: . dummy - optional user-defined context, as set by SNESSetFunction()
381: Output Parameter:
382: . F - function vector
383: */
384: static PetscErrorCode FormSource(TS ts, PetscReal ftime, Vec X_dummmy, Vec F, void *dummy)
385: {
386: PetscReal new_imp_rate;
387: LandauCtx *ctx;
388: DM pack;
389: REctx *rectx;
392: TSGetDM(ts,&pack);
393: DMGetApplicationContext(pack, &ctx);
394: rectx = (REctx*)ctx->data;
395: /* check for impurities */
396: rectx->impuritySrcRate(ftime,&new_imp_rate,ctx);
397: if (new_imp_rate != 0) {
398: if (new_imp_rate != rectx->current_rate) {
399: PetscInt ii;
400: PetscReal dne_dt,dni_dt,tilda_ns[LANDAU_MAX_SPECIES],temps[LANDAU_MAX_SPECIES];
401: Vec globFarray[LANDAU_MAX_GRIDS*LANDAU_MAX_BATCH_SZ];
402: rectx->current_rate = new_imp_rate;
403: for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) tilda_ns[ii] = 0;
404: for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) temps[ii] = 1;
405: dni_dt = new_imp_rate /* *ctx->t_0 */; /* fully ionized immediately, no normalize, stay in non-dim */
406: dne_dt = new_imp_rate*rectx->Ne_ion /* *ctx->t_0 */;
407: tilda_ns[0] = dne_dt; tilda_ns[rectx->imp_idx] = dni_dt;
408: temps[0] = rectx->T_cold; temps[rectx->imp_idx] = rectx->T_cold;
409: PetscInfo(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);
410: DMCompositeGetAccessArray(pack, F, ctx->num_grids*ctx->batch_sz, NULL, globFarray);
411: for (PetscInt grid=0 ; grid<ctx->num_grids ; grid++) {
412: /* add it */
413: DMPlexLandauAddMaxwellians(ctx->plex[grid],globFarray[ LAND_PACK_IDX(0,grid) ],ftime,temps,tilda_ns,grid,0,ctx);
414: VecViewFromOptions(globFarray[ LAND_PACK_IDX(0,grid) ],NULL,"-vec_view_sources");
415: }
416: // Does DMCompositeRestoreAccessArray copy the data back? (no)
417: DMCompositeRestoreAccessArray(pack, F, ctx->num_grids*ctx->batch_sz, NULL, globFarray);
418: }
419: } else {
420: VecZeroEntries(F);
421: rectx->current_rate = 0;
422: }
423: return 0;
424: }
425: PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx)
426: {
427: LandauCtx *ctx = (LandauCtx*) actx; /* user-defined application context */
428: REctx *rectx = (REctx*)ctx->data;
429: DM pack;
430: Vec globXArray[LANDAU_MAX_GRIDS*LANDAU_MAX_BATCH_SZ];
431: TSConvergedReason reason;
433: VecGetDM(X, &pack);
434: DMCompositeGetAccessArray(pack, X, ctx->num_grids*ctx->batch_sz, NULL, globXArray);
435: if (stepi > rectx->plotStep && rectx->plotting) {
436: rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */
437: rectx->plotIdx++;
438: }
439: /* view */
440: TSGetConvergedReason(ts,&reason);
441: if (time/rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) {
442: if ((reason || stepi==0 || rectx->plotIdx%rectx->print_period==0) && ctx->verbose > 0) {
443: /* print norms */
444: DMPlexLandauPrintNorms(X, stepi);
445: }
446: if (!rectx->plotting) { /* first step of possible backtracks */
447: rectx->plotting = PETSC_TRUE;
448: /* diagnostics + change E field with Sptizer (not just a monitor) */
449: rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
450: } else {
451: PetscPrintf(PETSC_COMM_WORLD, "\t\t ERROR SKIP test spit ------\n");
452: rectx->plotting = PETSC_TRUE;
453: }
454: PetscObjectSetName((PetscObject) globXArray[ LAND_PACK_IDX(ctx->batch_view_idx,rectx->grid_view_idx) ], rectx->grid_view_idx==0 ? "ue" : "ui");
455: /* view, overwrite step when back tracked */
456: DMSetOutputSequenceNumber(pack, rectx->plotIdx, time*ctx->t_0);
457: VecViewFromOptions(globXArray[ LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx) ],NULL,"-vec_view");
459: rectx->plotStep = stepi;
460: } else {
461: if (rectx->plotting) PetscPrintf(PETSC_COMM_WORLD," ERROR rectx->plotting=%D step %D\n",rectx->plotting,stepi);
462: /* diagnostics + change E field with Sptizer (not just a monitor) - can we lag this? */
463: rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
464: }
465: /* parallel check that only works of all batches are identical */
466: if (reason && ctx->verbose > 3) {
467: PetscReal val,rval;
468: PetscMPIInt rank;
469: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
470: for (PetscInt grid=0;grid<ctx->num_grids;grid++) {
471: PetscInt nerrors=0;
472: for (PetscInt i=0; i<ctx->batch_sz;i++) {
473: VecNorm(globXArray[ LAND_PACK_IDX(i,grid) ],NORM_2,&val);
474: if (i==0) rval = val;
475: else if ((val=PetscAbs(val-rval)/rval) > 1000*PETSC_MACHINE_EPSILON) {
476: PetscPrintf(PETSC_COMM_SELF, " [%D] Warning %D.%D) diff = %2.15e\n",rank,grid,i,val);
477: nerrors++;
478: }
479: }
480: if (nerrors) {
481: PetscPrintf(PETSC_COMM_SELF, " ***** [%D] ERROR max %D errors\n",rank,nerrors);
482: } else {
483: PetscPrintf(PETSC_COMM_WORLD, "[%D] %D) batch consistency check OK\n",rank,grid);
484: }
485: }
486: }
487: rectx->idx = 0;
488: DMCompositeRestoreAccessArray(pack, X, ctx->num_grids*ctx->batch_sz, NULL, globXArray);
489: return 0;
490: }
492: PetscErrorCode PreStep(TS ts)
493: {
494: LandauCtx *ctx;
495: REctx *rectx;
496: DM dm;
497: PetscInt stepi;
498: PetscReal time;
499: Vec X;
502: /* not used */
503: TSGetDM(ts,&dm);
504: TSGetTime(ts,&time);
505: TSGetSolution(ts,&X);
506: DMGetApplicationContext(dm, &ctx);
507: rectx = (REctx*)ctx->data;
508: TSGetStepNumber(ts, &stepi);
509: /* update E */
510: rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez);
511: return 0;
512: }
514: /* 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) */
515: static PetscErrorCode stepSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
516: {
517: REctx *rectx = (REctx*)ctx->data;
520: if (time >= rectx->pulse_start) *rho = rectx->pulse_rate;
521: else *rho = 0.;
522: return 0;
523: }
524: static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
525: {
527: *rho = 0.;
528: return 0;
529: }
530: static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
531: {
532: REctx *rectx = (REctx*)ctx->data;
536: if (time < rectx->pulse_start || time > rectx->pulse_start + 3*rectx->pulse_width) *rho = 0;
537: /* else if (0) { */
538: /* 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); */
539: /* *rho = rectx->pulse_rate * (cycle / (stop - start)) / (1 + PetscExpReal(steep*(PetscSinReal(2*PETSC_PI*((t - start)/cycle + xi)) - PetscSinReal(2*PETSC_PI*xi)))); */
540: /* } else if (0) { */
541: /* double x = 2*(time - rectx->pulse_start)/(3*rectx->pulse_width) - 1; */
542: /* if (x==1 || x==-1) *rho = 0; */
543: /* else *rho = rectx->pulse_rate * PetscExpReal(-1/(1-x*x)); */
544: /* } */
545: else {
546: double x = PetscSinReal((time-rectx->pulse_start)/(3*rectx->pulse_width)*2*PETSC_PI - PETSC_PI/2) + 1; /* 0:2, integrates to 1.0 */
547: *rho = rectx->pulse_rate * x / (3*rectx->pulse_width);
548: if (!rectx->use_spitzer_eta) rectx->use_spitzer_eta = PETSC_TRUE; /* use it next time */
549: }
550: return 0;
551: }
555: static PetscErrorCode ProcessREOptions(REctx *rectx, const LandauCtx *ctx, DM dm, const char prefix[])
556: {
557: PetscErrorCode ierr;
558: PetscFunctionList plist = NULL, testlist = NULL, elist = NULL;
559: char pname[256],testname[256],ename[256];
560: DM dm_dummy;
561: PetscBool Connor_E = PETSC_FALSE;
564: DMCreate(PETSC_COMM_WORLD,&dm_dummy);
565: rectx->Ne_ion = 1; /* number of electrons given up by impurity ion */
566: rectx->T_cold = .005; /* kev */
567: rectx->ion_potential = 15; /* ev */
568: rectx->L = 2;
569: rectx->X_0 = NULL;
570: rectx->imp_idx = ctx->num_species - 1; /* default ionized impurity as last one */
571: rectx->pulse_start = PETSC_MAX_REAL;
572: rectx->pulse_width = 1;
573: rectx->plotStep = PETSC_MAX_INT;
574: rectx->pulse_rate = 1.e-1;
575: rectx->current_rate = 0;
576: rectx->plotIdx = 0;
577: rectx->j = 0;
578: rectx->plotDt = 1.0;
579: rectx->plotting = PETSC_FALSE;
580: rectx->use_spitzer_eta = PETSC_FALSE;
581: rectx->idx = 0;
582: rectx->print_period = 10;
583: rectx->grid_view_idx = 0;
584: /* Register the available impurity sources */
585: PetscFunctionListAdd(&plist,"step",&stepSrc);
586: PetscFunctionListAdd(&plist,"none",&zeroSrc);
587: PetscFunctionListAdd(&plist,"pulse",&pulseSrc);
588: PetscStrcpy(pname,"none");
589: PetscFunctionListAdd(&testlist,"none",&testNone);
590: PetscFunctionListAdd(&testlist,"spitzer",&testSpitzer);
591: PetscFunctionListAdd(&testlist,"stable",&testStable);
592: PetscStrcpy(testname,"none");
593: PetscFunctionListAdd(&elist,"none",&ENone);
594: PetscFunctionListAdd(&elist,"induction",&EInduction);
595: PetscFunctionListAdd(&elist,"constant",&EConstant);
596: PetscStrcpy(ename,"constant");
598: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");
599: PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "ex2.c", rectx->plotDt, &rectx->plotDt, NULL);
600: if (rectx->plotDt < 0) rectx->plotDt = 1e30;
601: if (rectx->plotDt == 0) rectx->plotDt = 1e-30;
602: PetscOptionsInt("-ex2_print_period", "Plotting interval", "ex2.c", rectx->print_period, &rectx->print_period, NULL);
603: PetscOptionsInt("-ex2_grid_view_idx", "grid_view_idx", "ex2.c", rectx->grid_view_idx, &rectx->grid_view_idx, NULL);
605: PetscOptionsFList("-ex2_impurity_source_type","Name of impurity source to run","",plist,pname,pname,sizeof(pname),NULL);
606: PetscOptionsFList("-ex2_test_type","Name of test to run","",testlist,testname,testname,sizeof(testname),NULL);
607: PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL);
609: PetscOptionsFList("-ex2_e_field_type","Electric field type","",elist,ename,ename,sizeof(ename),NULL);
610: rectx->Ne_ion = -ctx->charges[rectx->imp_idx]/ctx->charges[0];
611: PetscOptionsReal("-ex2_t_cold","Temperature of cold electron and ions after ionization in keV","none",rectx->T_cold,&rectx->T_cold, NULL);
612: PetscOptionsReal("-ex2_pulse_start_time","Time at which pulse happens for 'pulse' source","none",rectx->pulse_start,&rectx->pulse_start, NULL);
613: PetscOptionsReal("-ex2_pulse_width_time","Width of pulse 'pulse' source","none",rectx->pulse_width,&rectx->pulse_width, NULL);
614: PetscOptionsReal("-ex2_pulse_rate","Number density of pulse for 'pulse' source","none",rectx->pulse_rate,&rectx->pulse_rate, NULL);
615: rectx->T_cold *= 1.16e7; /* convert to Kelvin */
616: PetscOptionsReal("-ex2_ion_potential","Potential to ionize impurity (should be array) in ev","none",rectx->ion_potential,&rectx->ion_potential, NULL);
617: PetscOptionsReal("-ex2_inductance","Inductance E feild","none",rectx->L,&rectx->L, NULL);
618: PetscOptionsBool("-ex2_connor_e_field_units","Scale Ex but Connor-Hastie E_c","none",Connor_E,&Connor_E, NULL);
619: 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",rectx->Ne_ion,rectx->T_cold,rectx->ion_potential,ctx->Ez,ctx->v_0);
620: PetscOptionsEnd();
621: /* get impurity source rate function */
622: PetscFunctionListFind(plist,pname,&rectx->impuritySrcRate);
624: PetscFunctionListFind(testlist,testname,&rectx->test);
626: PetscFunctionListFind(elist,ename,&rectx->E);
628: PetscFunctionListDestroy(&plist);
629: PetscFunctionListDestroy(&testlist);
630: PetscFunctionListDestroy(&elist);
632: /* convert E from Connor-Hastie E_c units to real if doing Spitzer E */
633: if (Connor_E) {
634: PetscReal E = ctx->Ez, Tev = ctx->thermal_temps[0]*8.621738e-5, n = ctx->n_0*ctx->n[0];
635: CalculateE(Tev, n, ctx->lnLam, ctx->epsilon0, &E);
636: ((LandauCtx*)ctx)->Ez *= E;
637: }
638: DMDestroy(&dm_dummy);
639: return 0;
640: }
642: int main(int argc, char **argv)
643: {
644: DM pack;
645: Vec X,*XsubArray;
646: PetscInt dim = 2, nDMs;
647: TS ts;
648: Mat J;
649: PetscDS prob;
650: LandauCtx *ctx;
651: REctx *rectx;
652: #if defined PETSC_USE_LOG
653: PetscLogStage stage;
654: #endif
655: PetscMPIInt rank;
656: #if defined(PETSC_HAVE_THREADSAFETY)
657: double starttime, endtime;
658: #endif
659: PetscInitialize(&argc, &argv, NULL,help);
660: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
661: if (rank) { /* turn off output stuff for duplicate runs */
662: PetscOptionsClearValue(NULL,"-dm_view");
663: PetscOptionsClearValue(NULL,"-vec_view");
664: PetscOptionsClearValue(NULL,"-dm_view_diff");
665: PetscOptionsClearValue(NULL,"-vec_view_diff");
666: PetscOptionsClearValue(NULL,"-dm_view_sources");
667: PetscOptionsClearValue(NULL,"-vec_view_0");
668: PetscOptionsClearValue(NULL,"-dm_view_0");
669: PetscOptionsClearValue(NULL,"-vec_view_sources");
670: PetscOptionsClearValue(NULL,"-info"); /* this does not work */
671: }
672: PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);
673: /* Create a mesh */
674: DMPlexLandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, "", &X, &J, &pack);
675: DMCompositeGetNumberDM(pack,&nDMs);
676: PetscMalloc(sizeof(*XsubArray)*nDMs, &XsubArray);
677: PetscObjectSetName((PetscObject)J, "Jacobian");
678: PetscObjectSetName((PetscObject)X, "f");
679: DMGetApplicationContext(pack, &ctx);
680: DMSetUp(pack);
681: /* context */
682: PetscNew(&rectx);
683: ctx->data = rectx;
684: ProcessREOptions(rectx,ctx,pack,"");
685: DMGetDS(pack, &prob);
686: DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
687: PetscObjectSetName((PetscObject) XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx) ], rectx->grid_view_idx==0 ? "ue" : "ui");
688: DMViewFromOptions(ctx->plex[rectx->grid_view_idx],NULL,"-dm_view");
689: DMViewFromOptions(ctx->plex[rectx->grid_view_idx], NULL,"-dm_view_0");
690: VecViewFromOptions(XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx,rectx->grid_view_idx) ], NULL,"-vec_view_0"); // initial condition (monitor plots after step)
691: DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray); // read only
692: PetscFree(XsubArray);
693: VecViewFromOptions(X, NULL,"-vec_view_global"); // initial condition (monitor plots after step)
694: DMSetOutputSequenceNumber(pack, 0, 0.0);
695: /* Create timestepping solver context */
696: TSCreate(PETSC_COMM_SELF,&ts);
697: TSSetDM(ts,pack);
698: TSSetIFunction(ts,NULL,DMPlexLandauIFunction,NULL);
699: TSSetIJacobian(ts,J,J,DMPlexLandauIJacobian,NULL);
700: TSSetRHSFunction(ts,NULL,FormSource,NULL);
701: TSSetFromOptions(ts);
702: TSSetSolution(ts,X);
703: TSSetApplicationContext(ts, ctx);
704: TSMonitorSet(ts,Monitor,ctx,NULL);
705: TSSetPreStep(ts,PreStep);
706: rectx->Ez_initial = ctx->Ez; /* cache for induction caclulation - applied E field */
707: if (1) { /* warm up an test just DMPlexLandauIJacobian */
708: Vec vec;
709: PetscInt nsteps;
710: PetscReal dt;
711: PetscLogStageRegister("Warmup", &stage);
712: PetscLogStagePush(stage);
713: VecDuplicate(X,&vec);
714: VecCopy(X,vec);
715: TSGetMaxSteps(ts,&nsteps);
716: TSGetTimeStep(ts,&dt);
717: TSSetMaxSteps(ts,1);
718: TSSolve(ts,X);
719: TSSetMaxSteps(ts,nsteps);
720: TSSetStepNumber(ts,0);
721: TSSetTime(ts,0);
722: TSSetTimeStep(ts,dt);
723: rectx->plotIdx = 0;
724: rectx->plotting = PETSC_FALSE;
725: PetscLogStagePop();
726: VecCopy(vec,X);
727: VecDestroy(&vec);
728: ctx->aux_bool = PETSC_FALSE; // flag for not a clean Jacobian
729: }
730: /* go */
731: PetscLogStageRegister("Solve", &stage);
732: ctx->stage = 0; // lets not use this stage
733: #if defined(PETSC_HAVE_THREADSAFETY)
734: ctx->stage = 1; // not set with thread safety
735: #endif
736: TSSetSolution(ts,X);
737: PetscLogStagePush(stage);
738: #if defined(PETSC_HAVE_THREADSAFETY)
739: starttime = MPI_Wtime();
740: #endif
741: TSSolve(ts,X);
742: PetscLogStagePop();
743: #if defined(PETSC_HAVE_THREADSAFETY)
744: endtime = MPI_Wtime();
745: ctx->times[LANDAU_EX2_TSSOLVE] += (endtime - starttime);
746: #endif
747: VecViewFromOptions(X, NULL,"-vec_view_global");
748: /* clean up */
749: DMPlexLandauDestroyVelocitySpace(&pack);
750: TSDestroy(&ts);
751: VecDestroy(&X);
752: PetscFree(rectx);
753: PetscFinalize();
754: return 0;
755: }
757: /*TEST
759: testset:
760: requires: p4est !complex double
761: output_file: output/ex2_0.out
762: 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 -ksp_type tfqmr -pc_type jacobi -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 0 -dm_landau_gpu_assembly true -dm_landau_batch_size 2
763: test:
764: suffix: cpu
765: args: -dm_landau_device_type cpu -ksp_pc_side right
766: test:
767: suffix: kokkos
768: requires: kokkos_kernels
769: args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos -ksp_pc_side right
770: test:
771: suffix: cuda
772: requires: cuda
773: args: -dm_landau_device_type cuda -dm_mat_type aijcusparse -dm_vec_type cuda -mat_cusparse_use_cpu_solve -ksp_pc_side right
774: test:
775: suffix: kokkos_batch
776: requires: kokkos_kernels
777: 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
778: test:
779: suffix: kokkos_batch_coo
780: requires: kokkos_kernels
781: 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 -dm_landau_coo_assembly
783: TEST*/