diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp index 486d1a15a09..f30de76e3aa 100644 --- a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp +++ b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp @@ -59,24 +59,24 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time, { amrex::ignore_unused(a_step); - // Fields have E^{n}, B^{n-1/2} + // Fields have Eg^{n}, Bg^{n-1/2} // Particles have p^{n} and x^{n}. // Save the values at the start of the time step, m_WarpX->SaveParticlesAtImplicitStepStart ( ); - // Save the fields at the start of the step - m_Eold.Copy( m_WarpX->getMultiLevelField(FieldType::Efield_fp), FieldType::Efield_fp ); - m_E.Copy(m_Eold); // initial guess for E + // Save electric field at the start of the step + m_Eold.Copy( FieldType::Efield_fp ); - // Compute Bfield at time n+1/2 + // Advance WarpX owned Bfield_fp to t_{n+1/2} m_WarpX->EvolveB(a_dt, DtType::Full); m_WarpX->ApplyMagneticFieldBCs(); const amrex::Real half_time = a_time + 0.5_rt*a_dt; - // Solve nonlinear system for E at t_{n+1/2} + // Solve nonlinear system for Eg at t_{n+1/2} // Particles will be advanced to t_{n+1/2} + m_E.Copy( m_Eold ); // initial guess for Eg m_nlsolver->Solve( m_E, m_Eold, half_time, a_dt ); // Update WarpX owned Efield_fp to t_{n+1/2} @@ -85,8 +85,8 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time, // Advance particles from time n+1/2 to time n+1 m_WarpX->FinishImplicitParticleUpdate(); - // Advance E fields from time n+1/2 to time n+1 - // Eg^{n+1} = 2.0*E_g^{n+1/2} - E_g^n + // Advance Eg from time n+1/2 to time n+1 + // Eg^{n+1} = 2.0*Eg^{n+1/2} - Eg^n m_E.linComb( 2._rt, m_E, -1._rt, m_Eold ); m_WarpX->SetElectricFieldAndApplyBCs( m_E ); @@ -107,6 +107,6 @@ void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS, // current state of the fields E and B. Deposit current density at time n+1/2. m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian ); - // RHS = cvac^2*0.5*dt*( curl(B^{n+1/2}) - mu0*J^{n+1/2} ) + // RHS = cvac^2*0.5*dt*( curl(Bg^{n+1/2}) - mu0*Jg^{n+1/2} ) m_WarpX->ImplicitComputeRHSE(0.5_rt*a_dt, a_RHS); } diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp index d2be0012478..07e3606ee8c 100644 --- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp +++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp @@ -80,15 +80,14 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time, { amrex::ignore_unused(a_step); - // Fields have E^{n} and B^{n} + // Fields have Eg^{n} and Bg^{n} // Particles have p^{n} and x^{n}. // Save the values at the start of the time step, m_WarpX->SaveParticlesAtImplicitStepStart ( ); - // Save the fields at the start of the step - m_Eold.Copy( m_WarpX->getMultiLevelField(FieldType::Efield_fp), FieldType::Efield_fp ); - m_E.Copy(m_Eold); // initial guess for E + // Save electric field at the start of the step + m_Eold.Copy( FieldType::Efield_fp ); const int num_levels = static_cast(m_Bold.size()); for (int lev = 0; lev < num_levels; ++lev) { @@ -101,8 +100,9 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time, const amrex::Real theta_time = a_time + m_theta*a_dt; - // Solve nonlinear system for E at t_{n+theta} + // Solve nonlinear system for Eg at t_{n+theta} // Particles will be advanced to t_{n+1/2} + m_E.Copy( m_Eold ); // initial guess for Eg m_nlsolver->Solve( m_E, m_Eold, theta_time, a_dt ); // Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta} @@ -111,7 +111,7 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time, // Advance particles from time n+1/2 to time n+1 m_WarpX->FinishImplicitParticleUpdate(); - // Advance E and B fields from time n+theta to time n+1 + // Advance Eg and Bg from time n+theta to time n+1 const amrex::Real new_time = a_time + a_dt; FinishFieldUpdate( new_time ); @@ -132,7 +132,7 @@ void ThetaImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS, // current state of the fields E and B. Deposit current density at time n+1/2. m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian ); - // RHS = cvac^2*m_theta*dt*( curl(B^{n+theta}) - mu0*J^{n+1/2} ) + // RHS = cvac^2*m_theta*dt*( curl(Bg^{n+theta}) - mu0*Jg^{n+1/2} ) m_WarpX->ImplicitComputeRHSE(m_theta*a_dt, a_RHS); } @@ -154,8 +154,8 @@ void ThetaImplicitEM::FinishFieldUpdate ( amrex::Real a_new_time ) { amrex::ignore_unused(a_new_time); - // Eg^{n+1} = (1/theta)*E_g^{n+theta} + (1-1/theta)*E_g^n - // Bg^{n+1} = (1/theta)*B_g^{n+theta} + (1-1/theta)*B_g^n + // Eg^{n+1} = (1/theta)*Eg^{n+theta} + (1-1/theta)*Eg^n + // Bg^{n+1} = (1/theta)*Bg^{n+theta} + (1-1/theta)*Bg^n const amrex::Real c0 = 1._rt/m_theta; const amrex::Real c1 = 1._rt - c0; diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H index 19515c7f100..7d9db577e00 100644 --- a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H +++ b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H @@ -75,24 +75,7 @@ public: [[nodiscard]] RT dotProduct( const WarpXSolverVec& a_X ) const; - inline - void Copy ( const amrex::Vector, 3 > >& a_field_vec, - const warpx::fields::FieldType a_field_type ) - { - - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - IsDefined(), - "WarpXSolverVec::Copy() called on undefined WarpXSolverVec"); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - a_field_type==m_field_type, - "WarpXSolverVec::Copy() called with vecs of different types"); - for (int lev=0; lev,3>>& solver_vec = a_solver_vec.getVec(); - Copy(solver_vec, a_solver_vec.getVecType()); + + const amrex::Vector,3>>& this_vec = a_solver_vec.getVec(); + for (int lev=0; lev,3>>& this_vec(a_WarpX->getMultiLevelField(a_field_type)); + const amrex::Vector,3>>& this_vec( + a_WarpX->getMultiLevelField(a_field_type)); m_field_vec.resize(m_num_amr_levels); - const int lev = 0; - for (int n=0; n<3; n++) { - const amrex::MultiFab& mf_model = *this_vec[lev][n]; - m_field_vec[lev][n] = std::make_unique( mf_model.boxArray(), mf_model.DistributionMap(), - mf_model.nComp(), amrex::IntVect::TheZeroVector() ); + for (int lev=0; lev( mf_model.boxArray(), + mf_model.DistributionMap(), + mf_model.nComp(), + amrex::IntVect::TheZeroVector() ); + } } m_field_type = a_field_type; m_is_defined = true; @@ -28,6 +32,25 @@ void WarpXSolverVec::Define ( WarpX* a_WarpX, SetDotMask(); } +void WarpXSolverVec::Copy ( const warpx::fields::FieldType a_field_type ) +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + IsDefined(), + "WarpXSolverVec::Copy() called on undefined WarpXSolverVec"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + a_field_type==m_field_type, + "WarpXSolverVec::Copy() called with vecs of different types"); + + const amrex::Vector,3>>& this_vec( + m_WarpX->getMultiLevelField(a_field_type)); + for (int lev=0; lev