Alamo
Expression.cpp
Go to the documentation of this file.
1#include "Expression.H"
2
3namespace BC
4{
5
6void
7Expression::FillBoundary (amrex::BaseFab<Set::Scalar> &a_in,
8 const amrex::Box &a_box,
9 int ngrow, int /*dcomp*/, int /*ncomp*/, Set::Scalar time,
10 Orientation face, const amrex::Mask * /*mask*/)
11{
12 const amrex::Real* DX = m_geom.CellSize();
13
14 Util::Assert(INFO,TEST(a_in.nComp() == (int)m_ncomp));
15
16 amrex::Box box = a_box;
17 box.grow(ngrow);
18 const amrex::Dim3 lo= amrex::lbound(m_geom.Domain()), hi = amrex::ubound(m_geom.Domain());
19
20 amrex::Array4<amrex::Real> const& in = a_in.array();
21
22 amrex::IndexType type = amrex::IndexType::TheCellType();
23
24 for (int n = 0; n < a_in.nComp(); n++)
25 amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k)
26 {
27 Set::Vector pos = Set::Position(i, j, k, m_geom, type);
28 Set::Scalar x = 0.0, y=0.0, z=0.0, t=time;
29 x = pos(0);
30 #if AMREX_SPACEDIM > 1
31 y = pos(1);
32 #if AMREX_SPACEDIM > 2
33 z = pos(2);
34 #endif
35 #endif
36
37
38 amrex::IntVect glevel;
39 AMREX_D_TERM(glevel[0] = std::max(std::min(0,i-lo.x),i-hi.x); ,
40 glevel[1] = std::max(std::min(0,j-lo.y),j-hi.y); ,
41 glevel[2] = std::max(std::min(0,k-lo.z),k-hi.z); );
42
43 if (glevel[0]<0 && (face == Orientation::xlo || face == Orientation::All)) // Left boundary
44 {
45 if (BCUtil::IsDirichlet(m_bc_type[Face::XLO][n]))
46 in(i,j,k,n) = m_bc_func[Face::XLO][n](x,y,z,t);
47 else if(BCUtil::IsNeumann(m_bc_type[Face::XLO][n]))
48 in(i,j,k,n) = in(i-glevel[0],j,k,n) - (m_bc_func[Face::XLO].size() > 0 ? m_bc_func[Face::XLO][n](x,y,z,t)*DX[0] : 0);
49 else if(BCUtil::IsReflectEven(m_bc_type[Face::XLO][n]))
50 in(i,j,k,n) = in(1-glevel[0],j,k,n);
51 else if(BCUtil::IsReflectOdd(m_bc_type[Face::XLO][n]))
52 in(i,j,k,n) = -in(1-glevel[0],j,k,n);
53 else if(BCUtil::IsPeriodic(m_bc_type[Face::XLO][n])) {}
54 else
55 Util::Abort(INFO, "Incorrect boundary conditions");
56 }
57 else if (glevel[0]>0 && (face == Orientation::xhi || face == Orientation::All)) // Right boundary
58 {
59 if (BCUtil::IsDirichlet(m_bc_type[Face::XHI][n]))
60 in(i,j,k,n) = m_bc_func[Face::XHI][n](x,y,z,t);
61 else if(BCUtil::IsNeumann(m_bc_type[Face::XHI][n]))
62 in(i,j,k,n) = in(i-glevel[0],j,k,n) - (m_bc_func[Face::XHI].size() > 0 ? m_bc_func[Face::XHI][n](x,y,z,t)*DX[0] : 0);
63 else if(BCUtil::IsReflectEven(m_bc_type[Face::XHI][n]))
64 in(i,j,k,n) = in(hi.x-glevel[0],j,k,n);
65 else if(BCUtil::IsReflectOdd(m_bc_type[Face::XHI][n]))
66 in(i,j,k,n) = -in(hi.x-glevel[0],j,k,n);
67 else if(BCUtil::IsPeriodic(m_bc_type[Face::XHI][n])) {}
68 else
69 Util::Abort(INFO, "Incorrect boundary conditions");
70 }
71
72 else if (glevel[1]<0 && (face == Orientation::ylo || face == Orientation::All)) // Bottom boundary
73 {
74 if (BCUtil::IsDirichlet(m_bc_type[Face::YLO][n]))
75 in(i,j,k,n) = m_bc_func[Face::YLO][n](x,y,z,t);
76 else if (BCUtil::IsNeumann(m_bc_type[Face::YLO][n]))
77 in(i,j,k,n) = in(i,j-glevel[1],k,n) - (m_bc_func[Face::YLO].size() > 0 ? m_bc_func[Face::YLO][n](x,y,z,t)*DX[1] : 0);
78 else if (BCUtil::IsReflectEven(m_bc_type[Face::YLO][n]))
79 in(i,j,k,n) = in(i,j-glevel[1],k,n);
80 else if (BCUtil::IsReflectOdd(m_bc_type[Face::YLO][n]))
81 in(i,j,k,n) = -in(i,j-glevel[1],k,n);
82 else if(BCUtil::IsPeriodic(m_bc_type[Face::YLO][n])) {}
83 else
84 Util::Abort(INFO, "Incorrect boundary conditions");
85 }
86 else if (glevel[1]>0 && (face == Orientation::yhi || face == Orientation::All)) // Top boundary
87 {
88 if (BCUtil::IsDirichlet(m_bc_type[Face::YHI][n]))
89 in(i,j,k,n) = m_bc_func[Face::YHI][n](x,y,z,t);
90 else if (BCUtil::IsNeumann(m_bc_type[Face::YHI][n]))
91 in(i,j,k,n) = in(i,j-glevel[1],k,n) - (m_bc_func[Face::YHI].size() > 0 ? m_bc_func[Face::YHI][n](x,y,z,t)*DX[1] : 0);
92 else if (BCUtil::IsReflectEven(m_bc_type[Face::YHI][n]))
93 in(i,j,k,n) = in(i,hi.y-glevel[1],k,n);
94 else if (BCUtil::IsReflectOdd(m_bc_type[Face::YHI][n]))
95 in(i,j,k,n) = -in(i,hi.y-glevel[1],k,n);
96 else if(BCUtil::IsPeriodic(m_bc_type[Face::YHI][n])) {}
97 else
98 Util::Abort(INFO, "Incorrect boundary conditions");
99 }
100
101#if AMREX_SPACEDIM>2
102 else if (glevel[2]<0 && (face == Orientation::zlo || face == Orientation::All))
103 {
104 if (BCUtil::IsDirichlet(m_bc_type[Face::ZLO][n]))
105 in(i,j,k,n) = m_bc_func[Face::ZLO][n](x,y,z,t);
106 else if (BCUtil::IsNeumann(m_bc_type[Face::ZLO][n]))
107 in(i,j,k,n) = in(i,j,k-glevel[2],n) - (m_bc_func[Face::ZLO].size() > 0 ? m_bc_func[Face::ZLO][n](x,y,z,t)*DX[2] : 0);
108 else if (BCUtil::IsReflectEven(m_bc_type[Face::ZLO][n]))
109 in(i,j,k,n) = in(i,j,1-glevel[2],n);
110 else if (BCUtil::IsReflectOdd(m_bc_type[Face::ZLO][n]))
111 in(i,j,k,n) = -in(i,j,1-glevel[2],n);
112 else if(BCUtil::IsPeriodic(m_bc_type[Face::ZLO][n])) {}
113 else Util::Abort(INFO, "Incorrect boundary conditions");
114 }
115 else if (glevel[2]>0 && (face == Orientation::zhi || face == Orientation::All))
116 {
117 if (BCUtil::IsDirichlet(m_bc_type[Face::ZHI][n]))
118 in(i,j,k,n) = m_bc_func[Face::ZHI][n](x,y,z,t);
119 else if(BCUtil::IsNeumann(m_bc_type[Face::ZHI][n]))
120 in(i,j,k,n) = in(i,j,k-glevel[2],n) - (m_bc_func[Face::ZHI].size() > 0 ? m_bc_func[Face::ZHI][n](x,y,z,t)*DX[2] : 0);
121 else if(BCUtil::IsReflectEven(m_bc_type[Face::ZHI][n]))
122 in(i,j,k,n) = in(i,j,hi.z-glevel[2],n);
123 else if(BCUtil::IsReflectOdd(m_bc_type[Face::ZHI][n]))
124 in(i,j,k,n) = -in(i,j,hi.z-glevel[2],n);
125 else if(BCUtil::IsPeriodic(m_bc_type[Face::ZHI][n])) {}
126 else Util::Abort(INFO, "Incorrect boundary conditions");
127 }
128#endif
129
130
131 });
132}
133
134amrex::BCRec
135Expression::GetBCRec()
136{
137 int bc_lo[BL_SPACEDIM] = {AMREX_D_DECL(m_bc_type[Face::XLO][0],m_bc_type[Face::YLO][0],m_bc_type[Face::XLO][0])};
138 int bc_hi[BL_SPACEDIM] = {AMREX_D_DECL(m_bc_type[Face::XHI][0],m_bc_type[Face::YHI][0],m_bc_type[Face::XHI][0])};
139
140 return amrex::BCRec(bc_lo,bc_hi);
141}
142
143amrex::Array<int,AMREX_SPACEDIM>
144Expression::IsPeriodic()
145{
146 return {AMREX_D_DECL(BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
147 BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
148 BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))};
149}
150amrex::Periodicity Expression::Periodicity () const
151{
152 return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(m_geom.Domain().length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
153 m_geom.Domain().length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
154 m_geom.Domain().length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
155}
156amrex::Periodicity Expression::Periodicity (const amrex::Box& b) {
157 return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(b.length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
158 b.length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
159 b.length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
160
161}
162
163
164}
t
std::time_t t
Definition FileNameParse.cpp:9
TEST
#define TEST(x)
Definition Util.H:22
INFO
#define INFO
Definition Util.H:21
BC::BC< Set::Scalar >::m_geom
amrex::Geometry m_geom
Definition BC.H:129
BC::Expression::m_bc_type
std::array< std::vector< int >, m_nfaces > m_bc_type
Definition Expression.H:95
BC::Expression::GetBCRec
amrex::BCRec GetBCRec() override
Definition Expression.cpp:135
BC::Expression::FillBoundary
virtual void FillBoundary(amrex::BaseFab< Set::Scalar > &in, const amrex::Box &box, int ngrow, int dcomp, int ncomp, amrex::Real time, Orientation face=Orientation::All, const amrex::Mask *mask=nullptr) override
Definition Expression.cpp:7
BC::Expression::IsPeriodic
virtual amrex::Array< int, AMREX_SPACEDIM > IsPeriodic() override
Definition Expression.cpp:144
BC::Expression::m_bc_func
std::array< std::vector< amrex::ParserExecutor< 4 > >, m_nfaces > m_bc_func
Definition Expression.H:97
BC::Expression::m_ncomp
unsigned int m_ncomp
Definition Expression.H:93
BC::Expression::Periodicity
virtual amrex::Periodicity Periodicity() const override
Definition Expression.cpp:150
BC::BCUtil::IsDirichlet
bool IsDirichlet(int bctype)
Definition BC.cpp:48
BC::BCUtil::IsPeriodic
bool IsPeriodic(int bctype)
Definition BC.cpp:33
BC::BCUtil::IsNeumann
bool IsNeumann(int bctype)
Definition BC.cpp:41
BC::BCUtil::IsReflectOdd
bool IsReflectOdd(int bctype)
Definition BC.cpp:59
BC::BCUtil::IsReflectEven
bool IsReflectEven(int bctype)
Definition BC.cpp:54
BC
Collection of boundary condition (BC) objects.
Definition BC.cpp:5
BC::Orientation
Orientation
Definition BC.H:31
BC::All
@ All
Definition BC.H:32
BC::AMREX_D_DECL
@ AMREX_D_DECL
Definition BC.H:33
Set::Scalar
amrex::Real Scalar
Definition Base.H:18
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:19
Set::Position
AMREX_FORCE_INLINE Vector Position(const int &i, const int &j, const int &k, const amrex::Geometry &geom, const amrex::IndexType &ixType)
Definition Base.H:120
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:225
Util::Assert
AMREX_FORCE_INLINE void Assert(std::string file, std::string func, int line, std::string smt, bool pass, Args const &... args)
Definition Util.H:55
Generated by doxygen 1.9.8