Line data Source code
1 : #ifndef OPERATOR_H
2 : #define OPERATOR_H
3 :
4 : #include <AMReX_MLNodeLinOp.H>
5 : #include <AMReX_MLCellLinOp.H>
6 : #include <AMReX_MultiFabUtil.H>
7 : #include <AMReX_BaseFab.H>
8 :
9 : #include "BC/BC.H"
10 :
11 : using namespace amrex;
12 :
13 : enum Grid { Cell, Node};
14 :
15 : /// \brief Documentation for operator namespace
16 : namespace Operator
17 : {
18 :
19 : template <Grid G> class Operator;
20 :
21 : //
22 : //
23 : // NODE-BASED OPERATOR
24 : //
25 : //
26 :
27 : template <>
28 : class Operator<Grid::Node> : public amrex::MLNodeLinOp
29 : {
30 : //
31 : // Public: Constructor/Destructor/Operators/defines
32 : //
33 : public :
34 749 : Operator () {}
35 : Operator (const amrex::Vector<amrex::Geometry>& a_geom,
36 : const amrex::Vector<amrex::BoxArray>& a_grids,
37 : const Vector<DistributionMapping>& a_dmap,
38 : const LPInfo& a_info = LPInfo(),
39 : const Vector<FabFactory<FArrayBox> const*>& a_factory = {});
40 : virtual ~Operator ();
41 : Operator (const Operator&) = delete;
42 : Operator (Operator&&) = delete;
43 : Operator& operator= (const Operator&) = delete;
44 : Operator& operator= (Operator&&) = delete;
45 : void define (const Vector<Geometry>& a_geom, const Vector<BoxArray>& a_grids,
46 : const Vector<DistributionMapping>& a_dmap,
47 : const LPInfo& a_info = LPInfo(),
48 : const Vector<FabFactory<FArrayBox> const*>& a_factory = {});
49 18744 : const Geometry& Geom (int amr_lev, int mglev=0) const noexcept { return m_geom[amr_lev][mglev]; }
50 :
51 : void Reflux(int crse_amrlev,
52 : MultiFab& res, const MultiFab& crse_sol, const MultiFab& crse_rhs,
53 : MultiFab& fine_res, MultiFab& fine_sol, const MultiFab& fine_rhs)
54 : {reflux(crse_amrlev, res, crse_sol, crse_rhs,fine_res, fine_sol, fine_rhs);}
55 : void Apply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const { Fapply(amrlev,mglev,out,in);}
56 :
57 0 : void SetOmega(Set::Scalar a_omega) {m_omega = a_omega;}
58 0 : virtual void SetAverageDownCoeffs(bool) {Util::Abort(INFO,"Not implemented!");}
59 0 : void SetNormalizeDDW(bool a_normalize_ddw) {m_normalize_ddw = a_normalize_ddw;}
60 :
61 : //
62 : // Public Utilty functions
63 : //
64 : public:
65 : void RegisterNewFab(amrex::Vector<amrex::MultiFab> &input);
66 : void RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> > &input);
67 : const amrex::FArrayBox & GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter &mfi) const;
68 0 : virtual void SetHomogeneous (bool) {};
69 : //
70 : // Pure Virtual: you MUST override these functions
71 : //
72 : protected:
73 : virtual void Fapply (int amrlev, int mglev,MultiFab& out,const MultiFab& in) const override =0;
74 : virtual void averageDownCoeffs () = 0;
75 : //
76 : // Virtual: you SHOULD override these functions
77 : //
78 : protected:
79 : virtual void Diagonal (bool recompute=false);
80 : virtual void Diagonal (int amrlev, int mglev, amrex::MultiFab& diag);
81 : //
82 : // Virtual: you CAN override these functions (but probably don't need to)
83 : //
84 : virtual void Fsmooth (int amrlev, int mglev, MultiFab& x,const MultiFab& b) const override;
85 : virtual void normalize (int amrlev, int mglev, MultiFab& mf) const override;
86 : virtual void reflux (int crse_amrlev, MultiFab& res, const MultiFab& crse_sol, const MultiFab& crse_rhs,
87 : MultiFab& fine_res, MultiFab& fine_sol, const MultiFab& fine_rhs) const override;
88 : //
89 : // Virtual: you CAN'T override these functions (they take care of AMReX business) and you SHOULDN'T
90 : //
91 : protected:
92 : virtual void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final;
93 : virtual void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const override final;
94 : virtual void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, const MultiFab& fine_sol, const MultiFab& fine_rhs) final;
95 : virtual void prepareForSolve () override;
96 16006 : virtual bool isSingular (int amrlev) const final {return (amrlev == 0) ? m_is_bottom_singular : false; }
97 29972 : virtual bool isBottomSingular () const final { return m_is_bottom_singular; }
98 : virtual void applyBC (int amrlev, int mglev, MultiFab& phi, BCMode bc_mode, amrex::MLLinOp::StateMode /**/, bool skip_fillboundary=false) const final;
99 : virtual void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final;
100 : public:
101 77192 : virtual int getNGrow(int /*alev*/=0,int /*mglev*/=0) const override final {return 2;}
102 : virtual void solutionResidual (int amrlev, MultiFab& resid, MultiFab& x, const MultiFab& b,
103 : const MultiFab* crse_bcdata=nullptr) override final;
104 : virtual void correctionResidual (int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b,
105 : BCMode bc_mode, const MultiFab* crse_bcdata=nullptr) override final;
106 :
107 :
108 : public:
109 : static void realFillBoundary(MultiFab &phi, const Geometry &geom);
110 :
111 : private:
112 : bool m_is_bottom_singular = false;
113 : bool m_masks_built = false;
114 : /// \todo we need to get rid of this
115 : // static constexpr amrex::IntVect AMREX_D_DECL(dx = {AMREX_D_DECL(1,0,0)},
116 : // dy = {AMREX_D_DECL(0,1,0)},
117 : // dz = {AMREX_D_DECL(0,0,1)});
118 : protected:
119 : int m_num_a_fabs = 0;
120 : bool m_diagonal_computed = false;
121 : amrex::Vector<amrex::Vector<amrex::Vector<amrex::MultiFab> > > m_a_coeffs;
122 : amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab> > > m_diag;
123 : Set::Scalar m_omega = 2./3.;
124 : bool m_normalize_ddw = false;
125 : };
126 :
127 :
128 :
129 : //
130 : //
131 : // CELL-BASED OPERATOR
132 : //
133 : //
134 :
135 :
136 : template<>
137 : class Operator<Grid::Cell> : public amrex::MLCellLinOp
138 : {
139 : public:
140 :
141 : friend class MLMG;
142 : friend class MLCGSolver;
143 :
144 : Operator ();
145 0 : virtual ~Operator () {};
146 :
147 : Operator (const Operator&) = delete;
148 : Operator (Operator&&) = delete;
149 : Operator& operator= (const Operator&) = delete;
150 : Operator& operator= (Operator&&) = delete;
151 :
152 : void define (amrex::Vector<amrex::Geometry> a_geom,
153 : const amrex::Vector<amrex::BoxArray>& a_grids,
154 : const amrex::Vector<amrex::DistributionMapping>& a_dmap,
155 : BC::BC<Set::Scalar>& a_bc,
156 : const amrex::LPInfo& a_info = amrex::LPInfo(),
157 : const amrex::Vector<amrex::FabFactory<amrex::FArrayBox> const*>& a_factory = {});
158 :
159 : // virtual void setLevelBC (int amrlev, const amrex::MultiFab* levelbcdata) override final;
160 :
161 : protected:
162 :
163 : BC::BC<Set::Scalar> *m_bc;
164 :
165 : amrex::Vector<std::unique_ptr<amrex::MLMGBndry> > m_bndry_sol;
166 : amrex::Vector<std::unique_ptr<amrex::BndryRegister> > m_crse_sol_br;
167 :
168 : amrex::Vector<std::unique_ptr<amrex::MLMGBndry> > m_bndry_cor;
169 : amrex::Vector<std::unique_ptr<amrex::BndryRegister> > m_crse_cor_br;
170 :
171 : // In case of agglomeration, coarse MG grids on amr level 0 are
172 : // not simply coarsened from fine MG grids. So we need to build
173 : // bcond and bcloc for each MG level.
174 : using RealTuple = std::array<amrex::Real,2*BL_SPACEDIM>;
175 : using BCTuple = std::array<amrex::BoundCond,2*BL_SPACEDIM>;
176 : class BndryCondLoc
177 : {
178 : public:
179 : BndryCondLoc (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);
180 : void setLOBndryConds (const amrex::Geometry& geom, const amrex::Real* dx,
181 : const amrex::Array<BCType,AMREX_SPACEDIM>& lobc,
182 : const amrex::Array<BCType,AMREX_SPACEDIM>& hibc,
183 : int ratio, const amrex::RealVect& a_loc);
184 : const BCTuple& bndryConds (const amrex::MFIter& mfi) const {
185 : return bcond[mfi];
186 : }
187 : const RealTuple& bndryLocs (const amrex::MFIter& mfi) const {
188 : return bcloc[mfi];
189 : }
190 : private:
191 : amrex::LayoutData<BCTuple> bcond;
192 : amrex::LayoutData<RealTuple> bcloc;
193 : };
194 : amrex::Vector<amrex::Vector<std::unique_ptr<BndryCondLoc> > > m_bcondloc;
195 :
196 : // used to save interpolation coefficients of the first interior cells
197 : mutable amrex::Vector<amrex::Vector<amrex::BndryRegister> > m_undrrelxr;
198 :
199 : // boundary cell flags for covered, not_covered, outside_domain
200 : amrex::Vector<amrex::Vector<std::array<amrex::MultiMask,2*AMREX_SPACEDIM> > > m_maskvals;
201 :
202 : //
203 : // functions
204 : //
205 :
206 : void updateSolBC (int amrlev, const amrex::MultiFab& crse_bcdata) const;
207 : void updateCorBC (int amrlev, const amrex::MultiFab& crse_bcdata) const;
208 :
209 : virtual void prepareForSolve () override;
210 :
211 :
212 : // PURE VIRTUAL METHODS
213 :
214 : virtual void Fapply (int amrlev, int mglev, amrex::MultiFab& out, const amrex::MultiFab& in) const override = 0;
215 : virtual void Fsmooth (int amrlev, int mglev, amrex::MultiFab& sol, const amrex::MultiFab& rsh, int redblack) const override = 0;
216 : virtual void FFlux (int amrlev, const MFIter& mfi,
217 : const Array<FArrayBox*,AMREX_SPACEDIM>& flux,
218 : const FArrayBox& sol, Location loc, const int face_only=0) const override = 0;
219 :
220 : void RegisterNewFab(amrex::Vector<amrex::MultiFab> &input);
221 : void RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> > &input);
222 : const amrex::FArrayBox & GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter &mfi) ;
223 :
224 0 : virtual bool isSingular (int /*amrlev*/) const final { return false; }
225 0 : virtual bool isBottomSingular () const final { return false; }
226 :
227 0 : virtual amrex::Real getAScalar () const final { return 0.0; }
228 0 : virtual amrex::Real getBScalar () const final { return 0.0; }
229 :
230 0 : virtual amrex::MultiFab const* getACoeffs (int /*amrlev*/, int /*mglev*/) const final { return nullptr;}
231 0 : virtual std::array<amrex::MultiFab const*,AMREX_SPACEDIM> getBCoeffs (int /*amrlev*/, int /*mglev*/) const final {
232 : std::array<amrex::MultiFab const*,AMREX_SPACEDIM> ret;
233 0 : AMREX_D_TERM(ret[0] = nullptr;, ret[1] = nullptr;,ret[2] = nullptr;);
234 0 : return ret;}
235 :
236 0 : virtual std::unique_ptr<amrex::MLLinOp> makeNLinOp (int /*grid_size*/) const final {
237 0 : Util::Abort("MLABecLaplacian::makeNLinOp: Not implmented");
238 0 : return std::unique_ptr<MLLinOp>{};
239 : }
240 :
241 : void averageDownCoeffs ();
242 : void averageDownCoeffsSameAmrLevel (amrex::Vector<amrex::MultiFab>&);
243 : const amrex::FArrayBox & GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter &mfi) const;
244 :
245 : static constexpr amrex::IntVect AMREX_D_DECL(dx = {AMREX_D_DECL(1,0,0)},
246 : dy = {AMREX_D_DECL(0,1,0)},
247 : dz = {AMREX_D_DECL(0,0,1)});
248 : private:
249 :
250 : int m_num_a_fabs = 0;
251 : amrex::Vector<amrex::Vector<amrex::Vector<amrex::MultiFab> > > m_a_coeffs;
252 :
253 : };
254 :
255 :
256 :
257 : }
258 :
259 :
260 : #endif
|
