Line data Source code
1 : #include "Expression.H"
2 :
3 : namespace BC
4 : {
5 :
6 : void
7 0 : 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 0 : const amrex::Real* DX = m_geom.CellSize();
13 :
14 0 : Util::Assert(INFO,TEST(a_in.nComp() == (int)m_ncomp));
15 :
16 0 : amrex::Box box = a_box;
17 0 : box.grow(ngrow);
18 0 : const amrex::Dim3 lo= amrex::lbound(m_geom.Domain()), hi = amrex::ubound(m_geom.Domain());
19 :
20 0 : amrex::Array4<amrex::Real> const& in = a_in.array();
21 :
22 0 : amrex::IndexType type = amrex::IndexType::TheCellType();
23 :
24 0 : for (int n = 0; n < a_in.nComp(); n++)
25 : {
26 0 : amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k)
27 : {
28 0 : Set::Vector pos = Set::Position(i, j, k, m_geom, type);
29 0 : Set::Scalar x = 0.0, y=0.0, z=0.0, t=time;
30 0 : x = pos(0);
31 : #if AMREX_SPACEDIM > 1
32 0 : y = pos(1);
33 : #if AMREX_SPACEDIM > 2
34 0 : z = pos(2);
35 : #endif
36 : #endif
37 :
38 :
39 0 : amrex::IntVect glevel;
40 0 : AMREX_D_TERM(glevel[0] = std::max(std::min(0,i-lo.x),i-hi.x); ,
41 : glevel[1] = std::max(std::min(0,j-lo.y),j-hi.y); ,
42 : glevel[2] = std::max(std::min(0,k-lo.z),k-hi.z); );
43 :
44 0 : if (glevel[0]<0 && (face == Orientation::xlo || face == Orientation::All)) // Left boundary
45 : {
46 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::XLO][n]))
47 0 : in(i,j,k,n) = m_bc_func[Face::XLO][n](x,y,z,t);
48 0 : else if(BCUtil::IsNeumann(m_bc_type[Face::XLO][n]))
49 0 : 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);
50 0 : else if(BCUtil::IsReflectEven(m_bc_type[Face::XLO][n]))
51 0 : in(i,j,k,n) = in(1-glevel[0],j,k,n);
52 0 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::XLO][n]))
53 0 : in(i,j,k,n) = -in(1-glevel[0],j,k,n);
54 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::XLO][n])) {}
55 : else
56 0 : Util::Abort(INFO, "Incorrect boundary conditions");
57 : }
58 0 : else if (glevel[0]>0 && (face == Orientation::xhi || face == Orientation::All)) // Right boundary
59 : {
60 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::XHI][n]))
61 0 : in(i,j,k,n) = m_bc_func[Face::XHI][n](x,y,z,t);
62 0 : else if(BCUtil::IsNeumann(m_bc_type[Face::XHI][n]))
63 0 : 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);
64 0 : else if(BCUtil::IsReflectEven(m_bc_type[Face::XHI][n]))
65 0 : in(i,j,k,n) = in(hi.x-glevel[0],j,k,n);
66 0 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::XHI][n]))
67 0 : in(i,j,k,n) = -in(hi.x-glevel[0],j,k,n);
68 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::XHI][n])) {}
69 : else
70 0 : Util::Abort(INFO, "Incorrect boundary conditions");
71 : }
72 0 : else if (glevel[1]<0 && (face == Orientation::ylo || face == Orientation::All)) // Bottom boundary
73 : {
74 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::YLO][n]))
75 0 : in(i,j,k,n) = m_bc_func[Face::YLO][n](x,y,z,t);
76 0 : else if (BCUtil::IsNeumann(m_bc_type[Face::YLO][n]))
77 0 : 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 0 : else if (BCUtil::IsReflectEven(m_bc_type[Face::YLO][n]))
79 0 : in(i,j,k,n) = in(i,j-glevel[1],k,n);
80 0 : else if (BCUtil::IsReflectOdd(m_bc_type[Face::YLO][n]))
81 0 : in(i,j,k,n) = -in(i,j-glevel[1],k,n);
82 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::YLO][n])) {}
83 : else
84 0 : Util::Abort(INFO, "Incorrect boundary conditions");
85 : }
86 0 : else if (glevel[1]>0 && (face == Orientation::yhi || face == Orientation::All)) // Top boundary
87 : {
88 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::YHI][n]))
89 0 : in(i,j,k,n) = m_bc_func[Face::YHI][n](x,y,z,t);
90 0 : else if (BCUtil::IsNeumann(m_bc_type[Face::YHI][n]))
91 0 : 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 0 : else if (BCUtil::IsReflectEven(m_bc_type[Face::YHI][n]))
93 0 : in(i,j,k,n) = in(i,hi.y-glevel[1],k,n);
94 0 : else if (BCUtil::IsReflectOdd(m_bc_type[Face::YHI][n]))
95 0 : in(i,j,k,n) = -in(i,hi.y-glevel[1],k,n);
96 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::YHI][n])) {}
97 : else
98 0 : Util::Abort(INFO, "Incorrect boundary conditions");
99 : }
100 : #if AMREX_SPACEDIM>2
101 0 : else if (glevel[2]<0 && (face == Orientation::zlo || face == Orientation::All))
102 : {
103 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::ZLO][n]))
104 0 : in(i,j,k,n) = m_bc_func[Face::ZLO][n](x,y,z,t);
105 0 : else if (BCUtil::IsNeumann(m_bc_type[Face::ZLO][n]))
106 0 : 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);
107 0 : else if (BCUtil::IsReflectEven(m_bc_type[Face::ZLO][n]))
108 0 : in(i,j,k,n) = in(i,j,1-glevel[2],n);
109 0 : else if (BCUtil::IsReflectOdd(m_bc_type[Face::ZLO][n]))
110 0 : in(i,j,k,n) = -in(i,j,1-glevel[2],n);
111 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::ZLO][n])) {}
112 0 : else Util::Abort(INFO, "Incorrect boundary conditions");
113 : }
114 0 : else if (glevel[2]>0 && (face == Orientation::zhi || face == Orientation::All))
115 : {
116 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::ZHI][n]))
117 0 : in(i,j,k,n) = m_bc_func[Face::ZHI][n](x,y,z,t);
118 0 : else if(BCUtil::IsNeumann(m_bc_type[Face::ZHI][n]))
119 0 : 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);
120 0 : else if(BCUtil::IsReflectEven(m_bc_type[Face::ZHI][n]))
121 0 : in(i,j,k,n) = in(i,j,hi.z-glevel[2],n);
122 0 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::ZHI][n]))
123 0 : in(i,j,k,n) = -in(i,j,hi.z-glevel[2],n);
124 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::ZHI][n])) {}
125 0 : else Util::Abort(INFO, "Incorrect boundary conditions");
126 : }
127 : #endif
128 :
129 0 : });
130 :
131 : // Average physical corner ghost cells from the face ghosts, matching
132 : // Constant BC behavior and avoiding undefined corner values in centered stencils.
133 0 : amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k)
134 : {
135 0 : 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) );
136 0 : 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) );
137 0 : 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) );
138 0 : 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) );
139 0 : });
140 : }
141 0 : }
142 :
143 : amrex::BCRec
144 0 : Expression::GetBCRec()
145 : {
146 0 : 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])};
147 0 : 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])};
148 :
149 0 : return amrex::BCRec(bc_lo,bc_hi);
150 : }
151 :
152 : amrex::Array<int,AMREX_SPACEDIM>
153 0 : Expression::IsPeriodic()
154 : {
155 0 : return {AMREX_D_DECL(BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
156 : BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
157 0 : BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))};
158 : }
159 0 : amrex::Periodicity Expression::Periodicity () const
160 : {
161 0 : return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(m_geom.Domain().length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
162 : m_geom.Domain().length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
163 0 : m_geom.Domain().length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
164 : }
165 0 : amrex::Periodicity Expression::Periodicity (const amrex::Box& b) {
166 0 : return amrex::Periodicity(amrex::IntVect(AMREX_D_DECL(b.length(0) * BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
167 : b.length(1) * BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
168 0 : b.length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
169 :
170 : }
171 :
172 :
173 : }
|
