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 "ImplicitSolver_fwd.H"
 
#include "FieldSolver/ImplicitSolvers/ImplicitOptions.H"
 
#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H"
#include "NonlinearSolvers/NonlinearSolverLibrary.H"
 
#include <AMReX_Array.H>
#include <AMReX_INT.H>
#include <AMReX_REAL.H>
#include <AMReX_LO_BCTYPES.H>
 
#include <map>
#include <string>

 
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, bool from_restart) = 0;
 
    [[nodiscard]] bool IsDefined () const { return m_is_defined; }
 
    virtual void PrintParameters () const = 0;
 
    void CreateParticleAttributes () const;
 
    bool DoParticleSuborbits () const { return m_particle_suborbits; }

    virtual int 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;

    void PreLinearSolve ();
 
    [[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;
    }

    int GetMassMatricesPCnComp (const int field_dir, const int space_dir) const
    {
        if      (field_dir == 0) { return m_ncomp_pc_xx[space_dir]; }
        else if (field_dir == 1) { return m_ncomp_pc_yy[space_dir]; }
        else if (field_dir == 2) { return m_ncomp_pc_zz[space_dir]; }
        else { return -1; }
    }
 
    // 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 ApplyMassMatrices (
        ablastr::fields::MultiLevelVectorField& a_out,
        const ablastr::fields::MultiLevelVectorField& a_in,
        const ablastr::fields::MultiLevelVectorField* a_in_ref = nullptr,
        const ablastr::fields::MultiLevelVectorField* a_baseline = nullptr,
        amrex::Real a_scale = 1.0,
        bool a_zero_out_first = false );
 
    void ComputeJfromMassMatrices (bool a_J_from_MM_only);

    void FinishMassMatrices ();

    void FinishImplicitParticleUpdate (amrex::Real time, int step);
 
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;
 
    int m_suborbit_warning_threshold = 5;
    int m_suborbit_statistics_interval = 100;
    std::map<std::string, amrex::Long> m_accumulated_suborbit_counts;
    int m_suborbit_statistics_start_step = -1;
 
    void AccumulateSuborbitStatistics (
        std::map<std::string, amrex::Long> const& local_suborbit_counts,
        int step);

    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;

    int m_mass_matrices_pc_width = -1;

    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::IntVect m_ncomp_pc_xx = amrex::IntVect{0};
    amrex::IntVect m_ncomp_pc_yy = amrex::IntVect{0};
    amrex::IntVect m_ncomp_pc_zz = amrex::IntVect{0};

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

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

    void InitializeMassMatrices ();
 
    void PrintBaseImplicitSolverParameters () const;

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

    void SyncMassMatricesPCAndApplyBCs ();

    void SetMassMatricesForPC ( amrex::Real a_theta_dt );

    void CumulateJ ();

    void SaveE ();

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