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