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