Alamo
Expression.cpp
Go to the documentation of this file.
1 #include "Expression.H"
2 
3 namespace BC
4 {
5 
6 void
7 Expression::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 
134 amrex::BCRec
135 Expression::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 
143 amrex::Array<int,AMREX_SPACEDIM>
144 Expression::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 }
150 amrex::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 }
156 amrex::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 }
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:121
BC::AMREX_D_DECL
@ AMREX_D_DECL
Definition: BC.H:34
BC::Expression::m_bc_type
std::array< std::vector< int >, m_nfaces > m_bc_type
Definition: Expression.H:88
TEST
#define TEST(x)
Definition: Util.H:21
t
std::time_t t
Definition: FileNameParse.cpp:12
BC::All
@ All
Definition: BC.H:33
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
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition: Base.H:20
BC::Expression::GetBCRec
amrex::BCRec GetBCRec() override
Definition: Expression.cpp:135
Set::Scalar
amrex::Real Scalar
Definition: Base.H:19
BC::BC< Set::Scalar >::m_geom
amrex::Geometry m_geom
Definition: BC.H:129
BC::BCUtil::IsDirichlet
bool IsDirichlet(int bctype)
Definition: BC.cpp:48
BC::Expression::Periodicity
virtual amrex::Periodicity Periodicity() const override
Definition: Expression.cpp:150
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:65
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:165
BC::BCUtil::IsNeumann
bool IsNeumann(int bctype)
Definition: BC.cpp:41
BC::Orientation
Orientation
Definition: BC.H:32
BC
Collection of boundary condition (BC) objects.
Definition: BC.cpp:4
BC::Expression::IsPeriodic
virtual amrex::Array< int, AMREX_SPACEDIM > IsPeriodic() override
Definition: Expression.cpp:144
BC::BCUtil::IsPeriodic
bool IsPeriodic(int bctype)
Definition: BC.cpp:33
Expression.H
BC::BCUtil::IsReflectEven
bool IsReflectEven(int bctype)
Definition: BC.cpp:54
INFO
#define INFO
Definition: Util.H:20
BC::Expression::m_ncomp
unsigned int m_ncomp
Definition: Expression.H:86
BC::BCUtil::IsReflectOdd
bool IsReflectOdd(int bctype)
Definition: BC.cpp:59
BC::Expression::m_bc_func
std::array< std::vector< amrex::ParserExecutor< 4 > >, m_nfaces > m_bc_func
Definition: Expression.H:90
Generated on Fri Aug 16 2024 16:30:29 for Alamo by   doxygen 1.8.17