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 0 : amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k)
26 : {
27 0 : Set::Vector pos = Set::Position(i, j, k, m_geom, type);
28 0 : Set::Scalar x = 0.0, y=0.0, z=0.0, t=time;
29 0 : x = pos(0);
30 : #if AMREX_SPACEDIM > 1
31 0 : y = pos(1);
32 : #if AMREX_SPACEDIM > 2
33 0 : z = pos(2);
34 : #endif
35 : #endif
36 :
37 :
38 0 : amrex::IntVect glevel;
39 0 : 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 0 : if (glevel[0]<0 && (face == Orientation::xlo || face == Orientation::All)) // Left boundary
44 : {
45 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::XLO][n]))
46 0 : in(i,j,k,n) = m_bc_func[Face::XLO][n](x,y,z,t);
47 0 : else if(BCUtil::IsNeumann(m_bc_type[Face::XLO][n]))
48 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);
49 0 : else if(BCUtil::IsReflectEven(m_bc_type[Face::XLO][n]))
50 0 : in(i,j,k,n) = in(1-glevel[0],j,k,n);
51 0 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::XLO][n]))
52 0 : in(i,j,k,n) = -in(1-glevel[0],j,k,n);
53 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::XLO][n])) {}
54 : else
55 0 : Util::Abort(INFO, "Incorrect boundary conditions");
56 : }
57 0 : else if (glevel[0]>0 && (face == Orientation::xhi || face == Orientation::All)) // Right boundary
58 : {
59 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::XHI][n]))
60 0 : in(i,j,k,n) = m_bc_func[Face::XHI][n](x,y,z,t);
61 0 : else if(BCUtil::IsNeumann(m_bc_type[Face::XHI][n]))
62 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);
63 0 : else if(BCUtil::IsReflectEven(m_bc_type[Face::XHI][n]))
64 0 : in(i,j,k,n) = in(hi.x-glevel[0],j,k,n);
65 0 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::XHI][n]))
66 0 : in(i,j,k,n) = -in(hi.x-glevel[0],j,k,n);
67 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::XHI][n])) {}
68 : else
69 0 : Util::Abort(INFO, "Incorrect boundary conditions");
70 : }
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 :
101 : #if AMREX_SPACEDIM>2
102 0 : else if (glevel[2]<0 && (face == Orientation::zlo || face == Orientation::All))
103 : {
104 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::ZLO][n]))
105 0 : in(i,j,k,n) = m_bc_func[Face::ZLO][n](x,y,z,t);
106 0 : else if (BCUtil::IsNeumann(m_bc_type[Face::ZLO][n]))
107 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);
108 0 : else if (BCUtil::IsReflectEven(m_bc_type[Face::ZLO][n]))
109 0 : in(i,j,k,n) = in(i,j,1-glevel[2],n);
110 0 : else if (BCUtil::IsReflectOdd(m_bc_type[Face::ZLO][n]))
111 0 : in(i,j,k,n) = -in(i,j,1-glevel[2],n);
112 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::ZLO][n])) {}
113 0 : else Util::Abort(INFO, "Incorrect boundary conditions");
114 : }
115 0 : else if (glevel[2]>0 && (face == Orientation::zhi || face == Orientation::All))
116 : {
117 0 : if (BCUtil::IsDirichlet(m_bc_type[Face::ZHI][n]))
118 0 : in(i,j,k,n) = m_bc_func[Face::ZHI][n](x,y,z,t);
119 0 : else if(BCUtil::IsNeumann(m_bc_type[Face::ZHI][n]))
120 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);
121 0 : else if(BCUtil::IsReflectEven(m_bc_type[Face::ZHI][n]))
122 0 : in(i,j,k,n) = in(i,j,hi.z-glevel[2],n);
123 0 : else if(BCUtil::IsReflectOdd(m_bc_type[Face::ZHI][n]))
124 0 : in(i,j,k,n) = -in(i,j,hi.z-glevel[2],n);
125 0 : else if(BCUtil::IsPeriodic(m_bc_type[Face::ZHI][n])) {}
126 0 : else Util::Abort(INFO, "Incorrect boundary conditions");
127 : }
128 : #endif
129 :
130 :
131 0 : });
132 0 : }
133 :
134 : amrex::BCRec
135 0 : Expression::GetBCRec()
136 : {
137 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])};
138 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])};
139 :
140 0 : return amrex::BCRec(bc_lo,bc_hi);
141 : }
142 :
143 : amrex::Array<int,AMREX_SPACEDIM>
144 0 : Expression::IsPeriodic()
145 : {
146 0 : return {AMREX_D_DECL(BCUtil::IsPeriodic(m_bc_type[Face::XLO][0]),
147 : BCUtil::IsPeriodic(m_bc_type[Face::YLO][0]),
148 0 : BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))};
149 : }
150 0 : amrex::Periodicity Expression::Periodicity () const
151 : {
152 0 : 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 0 : m_geom.Domain().length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
155 : }
156 0 : amrex::Periodicity Expression::Periodicity (const amrex::Box& b) {
157 0 : 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 0 : b.length(2) * BCUtil::IsPeriodic(m_bc_type[Face::ZLO][0]))));
160 :
161 : }
162 :
163 :
164 : }
|
