diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index b08b41b27df..cb2da6f34df 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -104,9 +104,9 @@ WarpX::Evolve (int numsteps) // on first step, push p by -0.5*dt for (int lev = 0; lev <= finest_level; ++lev) { - mypc->PushP(lev, -0.5*dt[lev], - *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], - *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); + //mypc->PushP(lev, -0.5*dt[lev], + // *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], + // *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } is_synchronized = false; } else { @@ -138,7 +138,6 @@ WarpX::Evolve (int numsteps) amrex::Print() << "Error: do_subcycling = " << do_subcycling << std::endl; amrex::Abort("Unsupported do_subcycling type"); } - if (num_mirrors>0){ applyMirrors(cur_time); // E : guard cells are NOT up-to-date @@ -162,7 +161,7 @@ WarpX::Evolve (int numsteps) } is_synchronized = true; } - +//////////////////////////////////////////////// #ifdef WARPX_USE_PY if (warpx_py_afterEsolve) warpx_py_afterEsolve(); #endif @@ -183,13 +182,13 @@ WarpX::Evolve (int numsteps) } myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - +///////////////////////////////////////////////////// bool move_j = is_synchronized || do_insitu; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. ShiftGalileanBoundary(); - +/////////////////////// int num_moved = MoveWindow(move_j); // Electrostatic solver: particles can move by an arbitrary number of cells diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp index e1136ea3dc8..0f66df0bdc8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp @@ -319,7 +319,14 @@ AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ const Complex CRhonew= Rhonew_arr(i,j,k); const Complex Jcoef = Jcoef_arr(i,j,k); - + fields(i,j,k,Idx::Ex_avg) = i; + fields(i,j,k,Idx::Ey_avg) = j; + fields(i,j,k,Idx::Ez_avg) = k; + fields(i,j,k,Idx::Bx_avg) = 0.; // 3.e-8*i; + fields(i,j,k,Idx::By_avg) = 0.; // 3.e-8*j; + fields(i,j,k,Idx::Bz_avg) = 0.; // 3.e-8*k; + +/* //Update E (see the original Galilean article) fields(i,j,k,Idx::Ex) = T2*C*Ex_old + T2*S_ck*c2*I*(ky*Bz_old - kz*By_old) @@ -363,7 +370,8 @@ AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ fields(i,j,k,Idx::Bz_avg) = Psi1*Bz_old + I*Psi2*(kx*Ey_old - ky*Ex_old) + A1*I*(kx*Jy - ky*Jx)*inv_ep0; - +*/ + // fields(i,j,k,Idx::Ex_avg) = fields(i,j,k,Idx::Ex); // fields(i,j,k,Idx::Ey_avg) = fields(i,j,k,Idx::Ey); // fields(i,j,k,Idx::Ez_avg) = fields(i,j,k,Idx::Ez); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 514bed3b72c..5786c7a000b 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -91,13 +91,29 @@ namespace { solver.BackwardTransform(*Bfield[2], Idx::Bz); if (solver.fft_do_time_averaging){ + /* solver.BackwardTransform(*Efield_avg[0], Idx::Ex_avg); solver.BackwardTransform(*Efield_avg[1], Idx::Ey_avg); solver.BackwardTransform(*Efield_avg[2], Idx::Ez_avg); + solver.BackwardTransform(*Bfield_avg[0], Idx::Bx_avg); + solver.BackwardTransform(*Bfield_avg[1], Idx::By_avg); + solver.BackwardTransform(*Bfield_avg[2], Idx::Bz_avg); + */ + + solver.BackwardTransform(*Efield[0], Idx::Ex_avg); + solver.BackwardTransform(*Efield[1], Idx::Ey_avg); + solver.BackwardTransform(*Efield[2], Idx::Ez_avg); + solver.BackwardTransform(*Bfield[0], Idx::Bx_avg); + solver.BackwardTransform(*Bfield[1], Idx::By_avg); + solver.BackwardTransform(*Bfield[2], Idx::Bz_avg); + solver.BackwardTransform(*Efield_avg[0], Idx::Ex_avg); + solver.BackwardTransform(*Efield_avg[1], Idx::Ey_avg); + solver.BackwardTransform(*Efield_avg[2], Idx::Ez_avg); solver.BackwardTransform(*Bfield_avg[0], Idx::Bx_avg); solver.BackwardTransform(*Bfield_avg[1], Idx::By_avg); solver.BackwardTransform(*Bfield_avg[2], Idx::Bz_avg); + } } } @@ -113,6 +129,51 @@ WarpX::PushPSATD (amrex::Real a_dt) if (do_pml && pml[lev]->ok()) { pml[lev]->PushPSATD(); } + MultiFab& Ex = *Efield_aux[lev][0]; + MultiFab& Ey = *Efield_aux[lev][0]; + MultiFab& Ez = *Efield_aux[lev][0]; + MultiFab& Ex_avg = *Efield_avg_aux[lev][0]; + MultiFab& Ey_avg = *Efield_avg_aux[lev][0]; + MultiFab& Ez_avg = *Efield_avg_aux[lev][0]; + MultiFab& Bx = *Bfield_aux[lev][0]; + MultiFab& By = *Bfield_aux[lev][0]; + MultiFab& Bz = *Bfield_aux[lev][0]; + MultiFab& Bx_avg = *Bfield_avg_aux[lev][0]; + MultiFab& By_avg = *Bfield_avg_aux[lev][0]; + MultiFab& Bz_avg = *Bfield_avg_aux[lev][0]; + DistributionMapping dm { Ex.boxArray(), ParallelDescriptor::NProcs() }; + MultiFab diff(Ex.boxArray(), dm, 1, Ex.nGrow()); + Print()<<"RIGHT AFTER PushPSATD\n"; + diff.setVal(0.,0,1,Ex.nGrow()); + diff.plus(Ex,0,1,Ex.nGrow()); + diff.minus(Ex_avg,0,1,Ex.nGrow()); + Print()<<"Ex.max(0) "< mf{Efield_fp[lev][0].get(),Efield_fp[lev][1].get(),Efield_fp[lev][2].get()}; amrex::FillBoundary(mf, period); + Vector mf_avg{Efield_avg_fp[lev][0].get(),Efield_avg_fp[lev][1].get(),Efield_avg_fp[lev][2].get()}; + amrex::FillBoundary(mf_avg, period); } else { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( ng <= Efield_fp[lev][0]->nGrowVect(), @@ -406,6 +408,9 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type, IntVect ng) Efield_fp[lev][0]->FillBoundary(ng, period); Efield_fp[lev][1]->FillBoundary(ng, period); Efield_fp[lev][2]->FillBoundary(ng, period); + Efield_avg_fp[lev][0]->FillBoundary(ng, period); + Efield_avg_fp[lev][1]->FillBoundary(ng, period); + Efield_avg_fp[lev][2]->FillBoundary(ng, period); } } else if (patch_type == PatchType::coarse) @@ -423,6 +428,8 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type, IntVect ng) if ( safe_guard_cells ) { Vector mf{Efield_cp[lev][0].get(),Efield_cp[lev][1].get(),Efield_cp[lev][2].get()}; amrex::FillBoundary(mf, cperiod); + Vector mf_avg{Efield_avg_cp[lev][0].get(),Efield_avg_cp[lev][1].get(),Efield_avg_cp[lev][2].get()}; + amrex::FillBoundary(mf_avg, cperiod); } else { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -431,6 +438,9 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type, IntVect ng) Efield_cp[lev][0]->FillBoundary(ng, cperiod); Efield_cp[lev][1]->FillBoundary(ng, cperiod); Efield_cp[lev][2]->FillBoundary(ng, cperiod); + Efield_avg_cp[lev][0]->FillBoundary(ng, cperiod); + Efield_avg_cp[lev][1]->FillBoundary(ng, cperiod); + Efield_avg_cp[lev][2]->FillBoundary(ng, cperiod); } } } @@ -460,6 +470,8 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type, IntVect ng) if ( safe_guard_cells ) { Vector mf{Bfield_fp[lev][0].get(),Bfield_fp[lev][1].get(),Bfield_fp[lev][2].get()}; amrex::FillBoundary(mf, period); + Vector mf_avg{Bfield_avg_fp[lev][0].get(),Bfield_avg_fp[lev][1].get(),Bfield_avg_fp[lev][2].get()}; + amrex::FillBoundary(mf_avg, period); } else { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( ng <= Bfield_fp[lev][0]->nGrowVect(), @@ -467,6 +479,9 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type, IntVect ng) Bfield_fp[lev][0]->FillBoundary(ng, period); Bfield_fp[lev][1]->FillBoundary(ng, period); Bfield_fp[lev][2]->FillBoundary(ng, period); + Bfield_avg_fp[lev][0]->FillBoundary(ng, period); + Bfield_avg_fp[lev][1]->FillBoundary(ng, period); + Bfield_avg_fp[lev][2]->FillBoundary(ng, period); } } else if (patch_type == PatchType::coarse) @@ -484,6 +499,8 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type, IntVect ng) if ( safe_guard_cells ){ Vector mf{Bfield_cp[lev][0].get(),Bfield_cp[lev][1].get(),Bfield_cp[lev][2].get()}; amrex::FillBoundary(mf, cperiod); + Vector mf_avg{Bfield_avg_cp[lev][0].get(),Bfield_avg_cp[lev][1].get(),Bfield_avg_cp[lev][2].get()}; + amrex::FillBoundary(mf_avg, cperiod); } else { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( ng <= Bfield_cp[lev][0]->nGrowVect(), @@ -491,6 +508,9 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type, IntVect ng) Bfield_cp[lev][0]->FillBoundary(ng, cperiod); Bfield_cp[lev][1]->FillBoundary(ng, cperiod); Bfield_cp[lev][2]->FillBoundary(ng, cperiod); + Bfield_avg_cp[lev][0]->FillBoundary(ng, cperiod); + Bfield_avg_cp[lev][1]->FillBoundary(ng, cperiod); + Bfield_avg_cp[lev][2]->FillBoundary(ng, cperiod); } } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 44fca164867..f18f5439293 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1015,18 +1015,53 @@ PhysicalParticleContainer::FieldGather (int lev, } void -PhysicalParticleContainer::Evolve (int lev, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, - const MultiFab& Ex_avg, const MultiFab& Ey_avg, const MultiFab& Ez_avg, - const MultiFab& Bx_avg, const MultiFab& By_avg, const MultiFab& Bz_avg, - MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, - MultiFab* rho, MultiFab* crho, - const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, - const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real /*t*/, Real dt, DtType a_dt_type) +PhysicalParticleContainer::Evolve ( + int lev, + const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, + const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + const MultiFab& Ex_avg, const MultiFab& Ey_avg, const MultiFab& Ez_avg, + const MultiFab& Bx_avg, const MultiFab& By_avg, const MultiFab& Bz_avg, + MultiFab& jx, MultiFab& jy, MultiFab& jz, + MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, + MultiFab* rho, MultiFab* crho, + const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, + const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, + Real /*t*/, Real dt, DtType a_dt_type) { + DistributionMapping dm { Ex.boxArray(), ParallelDescriptor::NProcs() }; + MultiFab diff(Ex.boxArray(), dm, 1, Ex.nGrow()); + Print()<<"BEFORE\n"; + diff.setVal(0.,0,1,Ex.nGrow()); + diff.plus(Ex,0,1,Ex.nGrow()); + diff.minus(Ex_avg,0,1,Ex.nGrow()); + Print()<<"Ex.max(0) "<