aboutsummaryrefslogtreecommitdiff
path: root/Source/WarpXEvolve.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/WarpXEvolve.cpp')
-rw-r--r--Source/WarpXEvolve.cpp169
1 files changed, 32 insertions, 137 deletions
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index 89a73a978..a92f033b5 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -13,137 +13,16 @@ void
WarpX::Evolve (int numsteps) {
BL_PROFILE_REGION("WarpX::Evolve()");
+#ifdef WARPX_DO_ELECTROSTATIC
if (do_electrostatic) {
EvolveES(numsteps);
} else {
EvolveEM(numsteps);
}
-}
-
-void
-WarpX::EvolveES (int numsteps) {
-
- amrex::Print() << "Running in electrostatic mode \n";
-
- BL_PROFILE("WarpX::EvolveES()");
- Real cur_time = t_new[0];
- static int last_plot_file_step = 0;
- static int last_check_file_step = 0;
-
- int numsteps_max;
- if (numsteps < 0) { // Note that the default argument is numsteps = -1
- numsteps_max = max_step;
- } else {
- numsteps_max = std::min(istep[0]+numsteps, max_step);
- }
-
- bool max_time_reached = false;
-
- // nodal storage for thee electrostatic case
- const int num_levels = max_level + 1;
- Vector<std::unique_ptr<MultiFab> > rhoNodal(num_levels);
- Vector<std::unique_ptr<MultiFab> > phiNodal(num_levels);
- Vector<std::array<std::unique_ptr<MultiFab>, 3> > eFieldNodal(num_levels);
- const int ng = 1;
- for (int lev = 0; lev <= max_level; lev++) {
- BoxArray nba = boxArray(lev);
- nba.surroundingNodes();
- rhoNodal[lev].reset(new MultiFab(nba, dmap[lev], 1, ng));
- phiNodal[lev].reset(new MultiFab(nba, dmap[lev], 1, 2));
-
- eFieldNodal[lev][0].reset(new MultiFab(nba, dmap[lev], 1, ng));
- eFieldNodal[lev][1].reset(new MultiFab(nba, dmap[lev], 1, ng));
- eFieldNodal[lev][2].reset(new MultiFab(nba, dmap[lev], 1, ng));
- }
-
- const int lev = 0;
- for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step)
- {
-
- // Start loop on time steps
- amrex::Print() << "\nSTEP " << step+1 << " starts ...\n";
-
- // At initialization, particles have p^{n-1/2} and x^{n-1/2}.
- // Beyond one step, particles have p^{n-1/2} and x^{n}.
- if (is_synchronized) {
- // on first step, push X by 0.5*dt
- mypc->PushXES(0.5*dt[lev]);
- UpdatePlasmaInjectionPosition(0.5*dt[lev]);
- mypc->Redistribute();
- mypc->DepositCharge(rhoNodal);
- computePhi(rhoNodal, phiNodal);
- computeE(eFieldNodal, phiNodal);
- is_synchronized = false;
- }
-
- mypc->FieldGatherES(eFieldNodal, gather_masks);
-
- const std::string& ppltfile = amrex::Concatenate("particles", istep[0], 5);
- auto& pc = mypc->GetParticleContainer(0);
- pc.WriteAsciiFile(ppltfile);
-
- // Evolve particles to p^{n+1/2} and x^{n+1}
- mypc->EvolveES(eFieldNodal, rhoNodal, cur_time, dt[lev]);
-
- mypc->DepositCharge(rhoNodal);
-
- computePhi(rhoNodal, phiNodal);
- computeE(eFieldNodal, phiNodal);
-
- if (cur_time + dt[0] >= stop_time - 1.e-3*dt[0] || step == numsteps_max-1) {
- // on last step, push by only 0.5*dt to synchronize all at n+1/2
- mypc->PushXES(-0.5*dt[lev]);
- UpdatePlasmaInjectionPosition(-0.5*dt[lev]);
- is_synchronized = true;
- }
-
- mypc->Redistribute();
-
- ++istep[0];
-
- cur_time += dt[0];
-
- bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0);
-
- amrex::Print()<< "STEP " << step+1 << " ends." << " TIME = "
- << cur_time << " DT = " << dt[0] << "\n";
-
- // sync up time
- for (int i = 0; i <= finest_level; ++i) {
- t_new[i] = cur_time;
- }
-
- if (to_make_plot) {
- // replace with ES field Gather
- mypc->DepositCharge(rhoNodal);
- computePhi(rhoNodal, phiNodal);
- phiNodal[0]->FillBoundary(Geom(0).periodicity());
- computeE(eFieldNodal, phiNodal);
- mypc->FieldGatherES(eFieldNodal, gather_masks);
- last_plot_file_step = step+1;
- WritePlotFileES(rhoNodal, phiNodal, eFieldNodal);
- }
-
- if (check_int > 0 && (step+1) % check_int == 0) {
- last_check_file_step = step+1;
- WriteCheckPointFile();
- }
-
- if (cur_time >= stop_time - 1.e-3*dt[0]) {
- max_time_reached = true;
- break;
- }
-
- // End loop on time steps
- }
-
- if (plot_int > 0 && istep[0] > last_plot_file_step && (max_time_reached || istep[0] >= max_step)) {
- WritePlotFileES(rhoNodal, phiNodal, eFieldNodal);
- }
-
- if (check_int > 0 && istep[0] > last_check_file_step && (max_time_reached || istep[0] >= max_step)) {
- WriteCheckPointFile();
- }
+#else
+ EvolveEM(numsteps);
+#endif // WARPX_DO_ELECTROSTATIC
+
}
void
@@ -402,7 +281,7 @@ WarpX::EvolveB (int lev, Real dt, DtType typ)
const Box& tbz = mfi.tilebox(Bz_nodal_flag);
// Call picsar routine for each tile
- WRPX_PXR_PUSH_BVEC(
+ warpx_push_bvec(
tbx.loVect(), tbx.hiVect(),
tby.loVect(), tby.hiVect(),
tbz.loVect(), tbz.hiVect(),
@@ -412,8 +291,8 @@ WarpX::EvolveB (int lev, Real dt, DtType typ)
BL_TO_FORTRAN_3D((*Bx)[mfi]),
BL_TO_FORTRAN_3D((*By)[mfi]),
BL_TO_FORTRAN_3D((*Bz)[mfi]),
- &dtsdx[0], &dtsdx[1], &dtsdx[2],
- &norder);
+ &dtsdx[0], &dtsdx[1], &dtsdx[2],
+ &WarpX::maxwell_fdtd_solver_id);
if (cost) {
Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
@@ -434,7 +313,9 @@ WarpX::EvolveB (int lev, Real dt, DtType typ)
const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp()
: pml[lev]->GetMultiSigmaBox_cp();
-
+ int patch_level = (ipatch == 0) ? lev : lev-1;
+ const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
+ const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
#ifdef _OPENMP
#pragma omp parallel
#endif
@@ -454,7 +335,9 @@ WarpX::EvolveB (int lev, Real dt, DtType typ)
BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
- WRPX_PML_SIGMA_STAR_TO_FORTRAN(sigba[mfi],dttype));
+ WRPX_PML_SIGMA_STAR_TO_FORTRAN(sigba[mfi],dttype),
+ &dtsdx[0], &dtsdx[1], &dtsdx[2],
+ &WarpX::maxwell_fdtd_solver_id);
}
}
}
@@ -531,7 +414,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
const Box& tez = mfi.tilebox(Ez_nodal_flag);
// Call picsar routine for each tile
- WRPX_PXR_PUSH_EVEC(
+ warpx_push_evec(
tex.loVect(), tex.hiVect(),
tey.loVect(), tey.hiVect(),
tez.loVect(), tez.hiVect(),
@@ -545,8 +428,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ)
BL_TO_FORTRAN_3D((*jy)[mfi]),
BL_TO_FORTRAN_3D((*jz)[mfi]),
&mu_c2_dt,
- &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
- &norder);
+ &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
if (F) {
WRPX_CLEAN_EVEC(tex.loVect(), tex.hiVect(),
@@ -729,9 +611,22 @@ void
WarpX::ComputeDt ()
{
const Real* dx = geom[max_level].CellSize();
- const Real deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]),
- + 1./(dx[1]*dx[1]),
- + 1./(dx[2]*dx[2]))) * PhysConst::c );
+ Real deltat = 0.;
+
+ if (maxwell_fdtd_solver_id == 0) {
+ // CFL time step Yee solver
+ deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]),
+ + 1./(dx[1]*dx[1]),
+ + 1./(dx[2]*dx[2]))) * PhysConst::c );
+ } else {
+ // CFL time step CKC solver
+#if (BL_SPACEDIM == 3)
+ const Real delta = std::min(dx[0],std::min(dx[1],dx[2]));
+#elif (BL_SPACEDIM == 2)
+ const Real delta = std::min(dx[0],dx[1]);
+#endif
+ deltat = cfl*delta/PhysConst::c;
+ }
dt.resize(0);
dt.resize(max_level+1,deltat);
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-21 18:48:36 +0000