Line data Source code
1 : #include "Constant.H"
2 :
3 : namespace BC
4 : {
5 :
6 :
7 0 : 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 0 : 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 0 : Util::Warning(INFO,"This method is going away. Please use pp.queryclass() instead.");
21 :
22 0 : m_ncomp = a_ncomp;
23 :
24 0 : m_bc_type[Face::XLO].resize(m_ncomp,BCUtil::ReadString(bc_lo_str[0]));
25 0 : m_bc_type[Face::XHI].resize(m_ncomp,BCUtil::ReadString(bc_hi_str[0]));
26 0 : m_bc_type[Face::YLO].resize(m_ncomp,BCUtil::ReadString(bc_lo_str[1]));
27 0 : m_bc_type[Face::YHI].resize(m_ncomp,BCUtil::ReadString(bc_hi_str[1]));
28 : #if AMREX_SPACEDIM == 3
29 0 : m_bc_type[Face::ZLO].resize(m_ncomp,BCUtil::ReadString(bc_lo_str[2]));
30 0 : m_bc_type[Face::ZHI].resize(m_ncomp,BCUtil::ReadString(bc_hi_str[2]));
31 : #endif
32 :
33 :
34 0 : m_bc_val[Face::XLO].resize(m_ncomp,NAN);
35 0 : m_bc_val[Face::XHI].resize(m_ncomp,NAN);
36 0 : m_bc_val[Face::YLO].resize(m_ncomp,NAN);
37 0 : m_bc_val[Face::YHI].resize(m_ncomp,NAN);
38 : #if AMREX_SPACEDIM == 3
39 0 : m_bc_val[Face::ZLO].resize(m_ncomp,NAN);
40 0 : m_bc_val[Face::ZHI].resize(m_ncomp,NAN);
41 : #endif
42 :
43 0 : for (unsigned int i=0;i<m_ncomp;i++)
44 : {
45 0 : if (_bc_lo_1.size() > 0) m_bc_val[Face::XLO][i] = _bc_lo_1[i];
46 0 : if (_bc_hi_1.size() > 0) m_bc_val[Face::XHI][i] = _bc_hi_1[i];
47 0 : if (_bc_lo_2.size() > 0) m_bc_val[Face::YLO][i] = _bc_lo_2[i];
48 0 : if (_bc_hi_2.size() > 0) m_bc_val[Face::YHI][i] = _bc_hi_2[i];
49 : #if AMREX_SPACEDIM == 3
50 0 : if (_bc_lo_3.size() > 0) m_bc_val[Face::ZLO][i] = _bc_lo_3[i];
51 0 : if (_bc_hi_3.size() > 0) m_bc_val[Face::ZHI][i] = _bc_hi_3[i];
52 : #endif
53 : }
54 0 : }
55 :
56 :
57 : //amrex::Mask& m
58 : void
59 914842 : 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 914842 : const amrex::Real* DX = m_geom.CellSize();
65 :
66 6403894 : Util::Assert(INFO,TEST(a_in.nComp() == (int)m_ncomp));
67 :
68 914842 : amrex::Box box = a_box;
69 914842 : box.grow(ngrow);
70 2744526 : const amrex::Dim3 lo= amrex::lbound(m_geom.Domain()), hi = amrex::ubound(m_geom.Domain());
71 :
72 914842 : amrex::Array4<amrex::Real> const& in = a_in.array();
73 :
74 2495810 : for (int n = 0; n < a_in.nComp(); n++)
75 : {
76 1580968 : amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k)
77 : {
78 855446214 : amrex::IntVect glevel;
79 2574812114 : AMREX_D_TERM( glevel[0] = std::max(std::min(0,i-lo.x),i-hi.x); ,
80 : glevel[1] = std::max(std::min(0,j-lo.y),j-hi.y); ,
81 : glevel[2] = std::max(std::min(0,k-lo.z),k-hi.z); );
82 :
83 855446214 : if (glevel[0]<0 && (face == Orientation::xlo || face == Orientation::All)) // Left boundary
84 : {
85 13101596 : if (BCUtil::IsDirichlet(m_bc_type[Face::XLO][n]))
86 18277472 : in(i,j,k,n) = m_bc_val[Face::XLO][n](time);
87 3962860 : else if(BCUtil::IsNeumann(m_bc_type[Face::XLO][n]))
88 14674000 : 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);
89 294360 : else if(BCUtil::IsReflectEven(m_bc_type[Face::XLO][n]))
90 0 : in(i,j,k,n) = in(1-glevel[0],j,k,n);
91 294360 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::XLO][n]))
92 0 : in(i,j,k,n) = -in(1-glevel[0],j,k,n);
93 294360 : else if(BCUtil::IsPeriodic(m_bc_type[Face::XLO][n])) {}
94 : else
95 0 : Util::Abort(INFO, "Incorrect boundary conditions");
96 : }
97 842344618 : else if (glevel[0]>0 && (face == Orientation::xhi || face == Orientation::All)) // Right boundary
98 : {
99 2203840 : if (BCUtil::IsDirichlet(m_bc_type[Face::XHI][n]))
100 2935760 : in(i,j,k,n) = m_bc_val[Face::XHI][n](time);
101 735960 : else if(BCUtil::IsNeumann(m_bc_type[Face::XHI][n]))
102 2492640 : 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);
103 112800 : else if(BCUtil::IsReflectEven(m_bc_type[Face::XHI][n]))
104 0 : in(i,j,k,n) = in(hi.x-glevel[0],j,k,n);
105 112800 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::XHI][n]))
106 0 : in(i,j,k,n) = -in(hi.x-glevel[0],j,k,n);
107 112800 : else if(BCUtil::IsPeriodic(m_bc_type[Face::XHI][n])) {}
108 : else
109 0 : Util::Abort(INFO, "Incorrect boundary conditions");
110 : }
111 :
112 840140778 : else if (glevel[1]<0 && (face == Orientation::ylo || face == Orientation::All)) // Bottom boundary
113 : {
114 29346942 : if (BCUtil::IsDirichlet(m_bc_type[Face::YLO][n]))
115 4344768 : in(i,j,k,n) = m_bc_val[Face::YLO][n](time);
116 27174558 : else if (BCUtil::IsNeumann(m_bc_type[Face::YLO][n]))
117 107624984 : 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);
118 268312 : else if (BCUtil::IsReflectEven(m_bc_type[Face::YLO][n]))
119 0 : in(i,j,k,n) = in(i,j-glevel[1],k,n);
120 268312 : else if (BCUtil::IsReflectOdd(m_bc_type[Face::YLO][n]))
121 0 : in(i,j,k,n) = -in(i,j-glevel[1],k,n);
122 268312 : else if(BCUtil::IsPeriodic(m_bc_type[Face::YLO][n])) {}
123 : else
124 0 : Util::Abort(INFO, "Incorrect boundary conditions");
125 : }
126 810793836 : else if (glevel[1]>0 && (face == Orientation::yhi || face == Orientation::All)) // Top boundary
127 : {
128 27563214 : if (BCUtil::IsDirichlet(m_bc_type[Face::YHI][n]))
129 2117168 : in(i,j,k,n) = m_bc_val[Face::YHI][n](time);
130 26504630 : else if (BCUtil::IsNeumann(m_bc_type[Face::YHI][n]))
131 105500792 : 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);
132 129432 : else if (BCUtil::IsReflectEven(m_bc_type[Face::YHI][n]))
133 0 : in(i,j,k,n) = in(i,hi.y-glevel[1],k,n);
134 129432 : else if (BCUtil::IsReflectOdd(m_bc_type[Face::YHI][n]))
135 0 : in(i,j,k,n) = -in(i,hi.y-glevel[1],k,n);
136 129432 : else if(BCUtil::IsPeriodic(m_bc_type[Face::YHI][n])) {}
137 : else
138 0 : Util::Abort(INFO, "Incorrect boundary conditions");
139 : }
140 :
141 : #if AMREX_SPACEDIM>2
142 6658048 : else if (glevel[2]<0 && (face == Orientation::zlo || face == Orientation::All))
143 : {
144 600256 : if (BCUtil::IsDirichlet(m_bc_type[Face::ZLO][n]))
145 1200512 : in(i,j,k,n) = m_bc_val[Face::ZLO][n](time);
146 0 : else if (BCUtil::IsNeumann(m_bc_type[Face::ZLO][n]))
147 0 : 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);
148 0 : else if (BCUtil::IsReflectEven(m_bc_type[Face::ZLO][n]))
149 0 : in(i,j,k,n) = in(i,j,1-glevel[2],n);
150 0 : else if (BCUtil::IsReflectOdd(m_bc_type[Face::ZLO][n]))
151 0 : in(i,j,k,n) = -in(i,j,1-glevel[2],n);
152 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::ZLO][n])) {}
153 0 : else Util::Abort(INFO, "Incorrect boundary conditions");
154 : }
155 6057792 : else if (glevel[2]>0 && (face == Orientation::zhi || face == Orientation::All))
156 : {
157 102144 : if (BCUtil::IsDirichlet(m_bc_type[Face::ZHI][n]))
158 204288 : in(i,j,k,n) = m_bc_val[Face::ZHI][n](time);
159 0 : else if(BCUtil::IsNeumann(m_bc_type[Face::ZHI][n]))
160 0 : 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);
161 0 : else if(BCUtil::IsReflectEven(m_bc_type[Face::ZHI][n]))
162 0 : in(i,j,k,n) = in(i,j,hi.z-glevel[2],n);
163 0 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::ZHI][n]))
164 0 : in(i,j,k,n) = -in(i,j,hi.z-glevel[2],n);
165 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::ZHI][n])) {}
166 0 : else Util::Abort(INFO, "Incorrect boundary conditions");
167 : }
168 : #endif
169 855446214 : });
170 :
171 : // Average the value of corner cells based on the neighboring ghost cells.
172 : // This fixes NAN issues that can arise from neumann conditions calculated in ghost cells
173 : // Fixes this issue in 2D only for now.
174 : //
175 : // TODO: more general fix for 3D cell-based fields
176 1580968 : amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k)
177 : {
178 856337412 : if (i < lo.x && j < lo.y) in(i,j,k,n) = 0.5*( in(i+1,j,k,n) + in(i,j+1,k,n) );
179 856191012 : if (i < lo.x && j > hi.y) in(i,j,k,n) = 0.5*( in(i+1,j,k,n) + in(i,j-1,k,n) );
180 855875538 : if (i > hi.x && j < lo.y) in(i,j,k,n) = 0.5*( in(i-1,j,k,n) + in(i,j+1,k,n) );
181 855644754 : if (i > hi.x && j > hi.y) in(i,j,k,n) = 0.5*( in(i-1,j,k,n) + in(i,j-1,k,n) );
182 855446214 : });
183 :
184 : }
185 914842 : }
186 :
187 : amrex::BCRec
188 61310 : Constant::GetBCRec()
189 : {
190 61310 : 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])};
191 61310 : 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])};
192 :
193 61310 : return amrex::BCRec(bc_lo,bc_hi);
194 : }
195 :
196 : amrex::Array<int,AMREX_SPACEDIM>
197 0 : Constant::IsPeriodic()
198 : {
199 0 : return {AMREX_D_DECL(BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
200 : BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
201 0 : BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))};
202 : }
203 0 : amrex::Periodicity Constant::Periodicity () const
204 : {
205 0 : return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(m_geom.Domain().length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
206 : m_geom.Domain().length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
207 0 : m_geom.Domain().length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
208 : }
209 0 : amrex::Periodicity Constant::Periodicity (const amrex::Box& b) {
210 0 : return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(b.length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
211 : b.length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
212 0 : b.length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
213 :
214 : }
215 :
216 :
217 : }
|
