diff options
Diffstat (limited to 'Source/Evolve/WarpXEvolveEM.cpp')
| -rw-r--r-- | Source/Evolve/WarpXEvolveEM.cpp | 27 |
1 files changed, 21 insertions, 6 deletions
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 32a4747db..c5dbd73a7 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -348,7 +348,7 @@ WarpX::OneStep_sub1 (Real curtime) RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); NodalSyncJ(fine_lev, PatchType::fine); - ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*nmodes); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -375,7 +375,7 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(coarse_lev, curtime); StoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, 1); + AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, nmodes); EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); @@ -402,7 +402,7 @@ WarpX::OneStep_sub1 (Real curtime) RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); NodalSyncJ(fine_lev, PatchType::fine); - ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*nmodes); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -428,7 +428,7 @@ WarpX::OneStep_sub1 (Real curtime) // by only half a coarse step (second half) RestoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, 1, 1); + AddRhoFromFineLevelandSumBoundary(coarse_lev, nmodes, nmodes); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); @@ -492,8 +492,23 @@ WarpX::ComputeDt () if (maxwell_fdtd_solver_id == 0) { // CFL time step Yee solver #ifdef WARPX_RZ - // Derived semi-analytically by R. Lehe - deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); + // In the rz case, the Courant limit has been evaluated + // semi-analytically by R. Lehe, and resulted in the following + // coefficients. For an explanation, see (not officially published) + // www.normalesup.org/~lehe/Disp_relation_Circ.pdf + // NB : Here the coefficient for m=1 as compared to this document, + // as it was observed in practice that this coefficient was not + // high enough (The simulation became unstable). + Real multimode_coeffs[6] = { 0.2105, 1.0, 3.5234, 8.5104, 15.5059, 24.5037 }; + Real multimode_alpha; + if (nmodes < 7) { + // Use the table of the coefficients + multimode_alpha = multimode_coeffs[nmodes-1]; + } else { + // Use a realistic extrapolation + multimode_alpha = (nmodes - 1)*(nmodes - 1) - 0.4; + } + deltat = cfl * 1./( std::sqrt((1+multimode_alpha)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); #else deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), + 1./(dx[1]*dx[1]), |