ImplicitSolver.H

Go to the documentation of this file.

/* Copyright 2024 Justin Angus, Debojyoti Ghosh
 *
 * This file is part of WarpX.
 *
 * License: BSD-3-Clause-LBNL
 */
#ifndef Implicit_Solver_H_
#define Implicit_Solver_H_
 
#include "FieldSolver/ImplicitSolvers/ImplicitOptions.H"
 
#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H"
#include "NonlinearSolvers/NonlinearSolverLibrary.H"
 
#include <AMReX_Array.H>
#include <AMReX_REAL.H>
#include <AMReX_LO_BCTYPES.H>

 
class WarpX;
enum class FieldBoundaryType;
class ImplicitSolver
{
public:
 
    ImplicitSolver() = default;
 
    virtual ~ImplicitSolver() = default;
 
    // Prohibit Move and Copy operations
    ImplicitSolver(const ImplicitSolver&) = delete;
    ImplicitSolver& operator=(const ImplicitSolver&) = delete;
    ImplicitSolver(ImplicitSolver&&) = delete;
    ImplicitSolver& operator=(ImplicitSolver&&) = delete;
 
    //
    // the following routines are called by WarpX
    //

    virtual void Define ( WarpX*  a_WarpX ) = 0;
 
    [[nodiscard]] bool IsDefined () const { return m_is_defined; }
 
    virtual void PrintParameters () const = 0;
 
    void CreateParticleAttributes () const;
 
    bool DoParticleSuborbits () { return m_particle_suborbits; }

    virtual void OneStep ( amrex::Real  a_time,
                           amrex::Real  a_dt,
                           int          a_step ) = 0;
 
    //
    // the following routines are called by the linear and nonlinear solvers
    //

    virtual void ComputeRHS ( WarpXSolverVec&  a_RHS,
                        const WarpXSolverVec&  a_E,
                              amrex::Real      a_time,
                              int              a_nl_iter,
                              bool             a_from_jacobian ) = 0;
 
    [[nodiscard]] int numAMRLevels () const { return m_num_amr_levels; }
 
    [[nodiscard]] const amrex::Geometry& GetGeometry (int) const;
    [[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryLo () const;
    [[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryHi () const;
    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCLo () const;
    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCHi () const;

    [[nodiscard]] const amrex::Vector<amrex::Array<amrex::MultiFab*,3>>* GetMassMatricesCoeff() const {
        if (m_use_mass_matrices_pc) { return &m_mmpc_mfarrvec; }
        else { return nullptr; }
    }

    [[nodiscard]] virtual const amrex::MultiFab* GetCurl2BCmask (const int lev, const int dir) const
    {
        amrex::ignore_unused(lev,dir);
        return nullptr;
    }
 
    // This parameter is used for the time-step fraction in the PC for implicit
    // treatment of light waves in the curl-curl MLMG solver.
    // This function should return zero if light waves are not treated implicitly
    [[nodiscard]] virtual amrex::Real  GetThetaForPC () const = 0;
 
    void ComputeJfromMassMatrices (const bool  a_J_from_MM_only);
 
protected:

    WarpX* m_WarpX;
 
    bool m_is_defined = false;

    int m_num_amr_levels = 1;

    mutable amrex::Real m_dt = 0.0;

    amrex::Real m_theta = 0.5;

    NonlinearSolverType m_nlsolver_type;
    std::unique_ptr<NonlinearSolver<WarpXSolverVec,ImplicitSolver>> m_nlsolver;

    amrex::ParticleReal m_particle_tolerance = amrex::ParticleReal(1.0e-10);

    int m_max_particle_iterations = 21;

    bool m_particle_suborbits = false;

    bool m_print_unconverged_particle_details = false;

    bool m_skip_particle_picard_init = false;

    bool m_use_mass_matrices = false;

    bool m_use_mass_matrices_jacobian = false;

    bool m_use_mass_matrices_pc = false;

    amrex::IntVect m_ncomp_xx = amrex::IntVect{0};
    amrex::IntVect m_ncomp_xy = amrex::IntVect{0};
    amrex::IntVect m_ncomp_xz = amrex::IntVect{0};
    amrex::IntVect m_ncomp_yx = amrex::IntVect{0};
    amrex::IntVect m_ncomp_yy = amrex::IntVect{0};
    amrex::IntVect m_ncomp_yz = amrex::IntVect{0};
    amrex::IntVect m_ncomp_zx = amrex::IntVect{0};
    amrex::IntVect m_ncomp_zy = amrex::IntVect{0};
    amrex::IntVect m_ncomp_zz = amrex::IntVect{0};

    amrex::Vector<amrex::Array<amrex::MultiFab*,3>> m_mmpc_mfarrvec;

    void parseNonlinearSolverParams ( const amrex::ParmParse&  pp );

    void InitializeMassMatrices ();
 
    void PrintBaseImplicitSolverParameters () const;

    void PreRHSOp ( const amrex::Real  a_cur_time,
                    const int          a_nl_iter,
                    const bool         a_from_jacobian );

    void SyncMassMatricesPCAndApplyBCs ();

    void SetMassMatricesForPC ( const amrex::Real a_theta_dt );

    void CumulateJ ();

    void SaveE ();

    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>&, const int bdry_side ) const;
 
};
 
#endif