Alamo
Constant.cpp
Go to the documentation of this file.
1 #include "Constant.H"
2 
3 namespace BC
4 {
5 
6 
7 Constant::Constant (int a_ncomp,
8  amrex::Vector<std::string> bc_hi_str,
9  amrex::Vector<std::string> bc_lo_str,
10  AMREX_D_DECL(amrex::Vector<amrex::Real> _bc_lo_1,
11  amrex::Vector<amrex::Real> _bc_lo_2,
12  amrex::Vector<amrex::Real> _bc_lo_3),
13  AMREX_D_DECL(amrex::Vector<amrex::Real> _bc_hi_1,
14  amrex::Vector<amrex::Real> _bc_hi_2,
15  amrex::Vector<amrex::Real> _bc_hi_3))
16  //:
17  //AMREX_D_DECL(bc_lo_1(_bc_lo_1),bc_lo_2(_bc_lo_2),bc_lo_3(_bc_lo_3)),
18  //AMREX_D_DECL(bc_hi_1(_bc_hi_1),bc_hi_2(_bc_hi_2),bc_hi_3(_bc_hi_3))
19 {
20  Util::Warning(INFO,"This method is going away. Please use pp.queryclass() instead.");
21 
22  m_ncomp = a_ncomp;
23 
24  m_bc_type[Face::XLO].resize(m_ncomp,BCUtil::ReadString(bc_lo_str[0]));
25  m_bc_type[Face::XHI].resize(m_ncomp,BCUtil::ReadString(bc_hi_str[0]));
26  m_bc_type[Face::YLO].resize(m_ncomp,BCUtil::ReadString(bc_lo_str[1]));
27  m_bc_type[Face::YHI].resize(m_ncomp,BCUtil::ReadString(bc_hi_str[1]));
28  #if AMREX_SPACEDIM == 3
29  m_bc_type[Face::ZLO].resize(m_ncomp,BCUtil::ReadString(bc_lo_str[2]));
30  m_bc_type[Face::ZHI].resize(m_ncomp,BCUtil::ReadString(bc_hi_str[2]));
31  #endif
32 
33 
34  m_bc_val[Face::XLO].resize(m_ncomp,NAN);
35  m_bc_val[Face::XHI].resize(m_ncomp,NAN);
36  m_bc_val[Face::YLO].resize(m_ncomp,NAN);
37  m_bc_val[Face::YHI].resize(m_ncomp,NAN);
38  #if AMREX_SPACEDIM == 3
39  m_bc_val[Face::ZLO].resize(m_ncomp,NAN);
40  m_bc_val[Face::ZHI].resize(m_ncomp,NAN);
41  #endif
42 
43  for (unsigned int i=0;i<m_ncomp;i++)
44  {
45  if (_bc_lo_1.size() > 0) m_bc_val[Face::XLO][i] = _bc_lo_1[i];
46  if (_bc_hi_1.size() > 0) m_bc_val[Face::XHI][i] = _bc_hi_1[i];
47  if (_bc_lo_2.size() > 0) m_bc_val[Face::YLO][i] = _bc_lo_2[i];
48  if (_bc_hi_2.size() > 0) m_bc_val[Face::YHI][i] = _bc_hi_2[i];
49  #if AMREX_SPACEDIM == 3
50  if (_bc_lo_3.size() > 0) m_bc_val[Face::ZLO][i] = _bc_lo_3[i];
51  if (_bc_hi_3.size() > 0) m_bc_val[Face::ZHI][i] = _bc_hi_3[i];
52  #endif
53  }
54 }
55 
56 
57 //amrex::Mask& m
58 void
59 Constant::FillBoundary (amrex::BaseFab<Set::Scalar> &a_in,
60  const amrex::Box &a_box,
61  int ngrow, int /*dcomp*/, int /*ncomp*/, amrex::Real time,
62  Orientation face, const amrex::Mask * /*mask*/)
63 {
64  const amrex::Real* DX = m_geom.CellSize();
65 
66  Util::Assert(INFO,TEST(a_in.nComp() == (int)m_ncomp));
67 
68  amrex::Box box = a_box;
69  box.grow(ngrow);
70  const amrex::Dim3 lo= amrex::lbound(m_geom.Domain()), hi = amrex::ubound(m_geom.Domain());
71 
72  amrex::Array4<amrex::Real> const& in = a_in.array();
73 
74  for (int n = 0; n < a_in.nComp(); n++)
75  amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k)
76  {
77  amrex::IntVect glevel;
78  AMREX_D_TERM(glevel[0] = std::max(std::min(0,i-lo.x),i-hi.x); ,
79  glevel[1] = std::max(std::min(0,j-lo.y),j-hi.y); ,
80  glevel[2] = std::max(std::min(0,k-lo.z),k-hi.z); );
81 
82  if (glevel[0]<0 && (face == Orientation::xlo || face == Orientation::All)) // Left boundary
83  {
84  if (BCUtil::IsDirichlet(m_bc_type[Face::XLO][n]))
85  in(i,j,k,n) = m_bc_val[Face::XLO][n](time);
86  else if(BCUtil::IsNeumann(m_bc_type[Face::XLO][n]))
87  in(i,j,k,n) = in(i-glevel[0],j,k,n) - (m_bc_val[Face::XLO].size() > 0 ? m_bc_val[Face::XLO][n](time)*DX[0] : 0);
88  else if(BCUtil::IsReflectEven(m_bc_type[Face::XLO][n]))
89  in(i,j,k,n) = in(1-glevel[0],j,k,n);
90  else if(BCUtil::IsReflectOdd(m_bc_type[Face::XLO][n]))
91  in(i,j,k,n) = -in(1-glevel[0],j,k,n);
92  else if(BCUtil::IsPeriodic(m_bc_type[Face::XLO][n])) {}
93  else
94  Util::Abort(INFO, "Incorrect boundary conditions");
95  }
96  else if (glevel[0]>0 && (face == Orientation::xhi || face == Orientation::All)) // Right boundary
97  {
98  if (BCUtil::IsDirichlet(m_bc_type[Face::XHI][n]))
99  in(i,j,k,n) = m_bc_val[Face::XHI][n](time);
100  else if(BCUtil::IsNeumann(m_bc_type[Face::XHI][n]))
101  in(i,j,k,n) = in(i-glevel[0],j,k,n) - (m_bc_val[Face::XHI].size() > 0 ? m_bc_val[Face::XHI][n](time)*DX[0] : 0);
102  else if(BCUtil::IsReflectEven(m_bc_type[Face::XHI][n]))
103  in(i,j,k,n) = in(hi.x-glevel[0],j,k,n);
104  else if(BCUtil::IsReflectOdd(m_bc_type[Face::XHI][n]))
105  in(i,j,k,n) = -in(hi.x-glevel[0],j,k,n);
106  else if(BCUtil::IsPeriodic(m_bc_type[Face::XHI][n])) {}
107  else
108  Util::Abort(INFO, "Incorrect boundary conditions");
109  }
110 
111  else if (glevel[1]<0 && (face == Orientation::ylo || face == Orientation::All)) // Bottom boundary
112  {
113  if (BCUtil::IsDirichlet(m_bc_type[Face::YLO][n]))
114  in(i,j,k,n) = m_bc_val[Face::YLO][n](time);
115  else if (BCUtil::IsNeumann(m_bc_type[Face::YLO][n]))
116  in(i,j,k,n) = in(i,j-glevel[1],k,n) - (m_bc_val[Face::YLO].size() > 0 ? m_bc_val[Face::YLO][n](time)*DX[1] : 0);
117  else if (BCUtil::IsReflectEven(m_bc_type[Face::YLO][n]))
118  in(i,j,k,n) = in(i,j-glevel[1],k,n);
119  else if (BCUtil::IsReflectOdd(m_bc_type[Face::YLO][n]))
120  in(i,j,k,n) = -in(i,j-glevel[1],k,n);
121  else if(BCUtil::IsPeriodic(m_bc_type[Face::YLO][n])) {}
122  else
123  Util::Abort(INFO, "Incorrect boundary conditions");
124  }
125  else if (glevel[1]>0 && (face == Orientation::yhi || face == Orientation::All)) // Top boundary
126  {
127  if (BCUtil::IsDirichlet(m_bc_type[Face::YHI][n]))
128  in(i,j,k,n) = m_bc_val[Face::YHI][n](time);
129  else if (BCUtil::IsNeumann(m_bc_type[Face::YHI][n]))
130  in(i,j,k,n) = in(i,j-glevel[1],k,n) - (m_bc_val[Face::YHI].size() > 0 ? m_bc_val[Face::YHI][n](time)*DX[1] : 0);
131  else if (BCUtil::IsReflectEven(m_bc_type[Face::YHI][n]))
132  in(i,j,k,n) = in(i,hi.y-glevel[1],k,n);
133  else if (BCUtil::IsReflectOdd(m_bc_type[Face::YHI][n]))
134  in(i,j,k,n) = -in(i,hi.y-glevel[1],k,n);
135  else if(BCUtil::IsPeriodic(m_bc_type[Face::YHI][n])) {}
136  else
137  Util::Abort(INFO, "Incorrect boundary conditions");
138  }
139 
140 #if AMREX_SPACEDIM>2
141  else if (glevel[2]<0 && (face == Orientation::zlo || face == Orientation::All))
142  {
143  if (BCUtil::IsDirichlet(m_bc_type[Face::ZLO][n]))
144  in(i,j,k,n) = m_bc_val[Face::ZLO][n](time);
145  else if (BCUtil::IsNeumann(m_bc_type[Face::ZLO][n]))
146  in(i,j,k,n) = in(i,j,k-glevel[2],n) - (m_bc_val[Face::ZLO].size() > 0 ? m_bc_val[Face::ZLO][n](time)*DX[2] : 0);
147  else if (BCUtil::IsReflectEven(m_bc_type[Face::ZLO][n]))
148  in(i,j,k,n) = in(i,j,1-glevel[2],n);
149  else if (BCUtil::IsReflectOdd(m_bc_type[Face::ZLO][n]))
150  in(i,j,k,n) = -in(i,j,1-glevel[2],n);
151  else if(BCUtil::IsPeriodic(m_bc_type[Face::ZLO][n])) {}
152  else Util::Abort(INFO, "Incorrect boundary conditions");
153  }
154  else if (glevel[2]>0 && (face == Orientation::zhi || face == Orientation::All))
155  {
156  if (BCUtil::IsDirichlet(m_bc_type[Face::ZHI][n]))
157  in(i,j,k,n) = m_bc_val[Face::ZHI][n](time);
158  else if(BCUtil::IsNeumann(m_bc_type[Face::ZHI][n]))
159  in(i,j,k,n) = in(i,j,k-glevel[2],n) - (m_bc_val[Face::ZHI].size() > 0 ? m_bc_val[Face::ZHI][n](time)*DX[2] : 0);
160  else if(BCUtil::IsReflectEven(m_bc_type[Face::ZHI][n]))
161  in(i,j,k,n) = in(i,j,hi.z-glevel[2],n);
162  else if(BCUtil::IsReflectOdd(m_bc_type[Face::ZHI][n]))
163  in(i,j,k,n) = -in(i,j,hi.z-glevel[2],n);
164  else if(BCUtil::IsPeriodic(m_bc_type[Face::ZHI][n])) {}
165  else Util::Abort(INFO, "Incorrect boundary conditions");
166  }
167 #endif
168 
169 
170  });
171 }
172 
173 amrex::BCRec
175 {
176  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])};
177  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])};
178 
179  return amrex::BCRec(bc_lo,bc_hi);
180 }
181 
182 amrex::Array<int,AMREX_SPACEDIM>
184 {
185  return {AMREX_D_DECL(BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
186  BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
187  BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))};
188 }
189 amrex::Periodicity Constant::Periodicity () const
190 {
191  return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(m_geom.Domain().length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
192  m_geom.Domain().length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
193  m_geom.Domain().length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
194 }
195 amrex::Periodicity Constant::Periodicity (const amrex::Box& b) {
196  return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(b.length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
197  b.length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
198  b.length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
199 
200 }
201 
202 
203 }
BC::Constant::m_bc_val
std::array< std::vector< Numeric::Interpolator::Linear< Set::Scalar > >, m_nfaces > m_bc_val
Definition: Constant.H:160
BC::AMREX_D_DECL
@ AMREX_D_DECL
Definition: BC.H:34
TEST
#define TEST(x)
Definition: Util.H:21
BC::All
@ All
Definition: BC.H:33
Constant.H
BC::Constant::m_bc_type
std::array< std::vector< int >, m_nfaces > m_bc_type
Definition: Constant.H:159
BC::Constant::IsPeriodic
virtual amrex::Array< int, AMREX_SPACEDIM > IsPeriodic() override
Definition: Constant.cpp:183
BC::Constant::m_ncomp
unsigned int m_ncomp
Definition: Constant.H:152
BC::Constant::Constant
Constant(int a_ncomp)
Definition: Constant.H:107
BC::Constant::GetBCRec
amrex::BCRec GetBCRec() override
Definition: Constant.cpp:174
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
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
BC::BCUtil::ReadString
int ReadString(std::string bcstring)
Definition: BC.cpp:8
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:165
BC::Constant::Periodicity
virtual amrex::Periodicity Periodicity() const override
Definition: Constant.cpp:189
BC::BCUtil::IsNeumann
bool IsNeumann(int bctype)
Definition: BC.cpp:41
BC::Constant::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: Constant.cpp:59
BC::Orientation
Orientation
Definition: BC.H:32
BC
Collection of boundary condition (BC) objects.
Definition: BC.cpp:4
BC::BCUtil::IsPeriodic
bool IsPeriodic(int bctype)
Definition: BC.cpp:33
Util::Warning
void Warning(std::string file, std::string func, int line, Args const &... args)
Definition: Util.H:171
BC::BCUtil::IsReflectEven
bool IsReflectEven(int bctype)
Definition: BC.cpp:54
INFO
#define INFO
Definition: Util.H:20
BC::BCUtil::IsReflectOdd
bool IsReflectOdd(int bctype)
Definition: BC.cpp:59