Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLABecLap_1D_K.H
Go to the documentation of this file.
1#ifndef AMREX_MLABECLAP_1D_K_H_
2#define AMREX_MLABECLAP_1D_K_H_
3#include <AMReX_Config.H>
4
5namespace amrex {
6
7template <typename T>
9void mlabeclap_adotx (int i, int, int, int n, Array4<T> const& y,
10 Array4<T const> const& x,
11 Array4<T const> const& a,
12 Array4<T const> const& bX,
13 GpuArray<T,AMREX_SPACEDIM> const& dxinv,
14 T alpha, T beta) noexcept
15{
16 const T dhx = beta*dxinv[0]*dxinv[0];
17 y(i,0,0,n) = alpha*a(i,0,0)*x(i,0,0,n)
18 - dhx * (bX(i+1,0,0,n)*(x(i+1,0,0,n) - x(i ,0,0,n))
19 - bX(i ,0,0,n)*(x(i ,0,0,n) - x(i-1,0,0,n)));
20}
21
22template <typename T>
24void mlabeclap_adotx_os (int i, int, int, int n, Array4<T> const& y,
25 Array4<T const> const& x,
26 Array4<T const> const& a,
27 Array4<T const> const& bX,
28 Array4<int const> const& osm,
29 GpuArray<T,AMREX_SPACEDIM> const& dxinv,
30 T alpha, T beta) noexcept
31{
32 const T dhx = beta*dxinv[0]*dxinv[0];
33 if (osm(i,0,0) == 0) {
34 y(i,0,0,n) = T(0.0);
35 } else {
36 y(i,0,0,n) = alpha*a(i,0,0)*x(i,0,0,n)
37 - dhx * (bX(i+1,0,0,n)*(x(i+1,0,0,n) - x(i ,0,0,n))
38 - bX(i ,0,0,n)*(x(i ,0,0,n) - x(i-1,0,0,n)));
39 }
40}
41
42template <typename T>
44void mlabeclap_normalize (int i, int, int, int n, Array4<T> const& x,
45 Array4<T const> const& a,
46 Array4<T const> const& bX,
47 GpuArray<T,AMREX_SPACEDIM> const& dxinv,
48 T alpha, T beta) noexcept
49{
50 const T dhx = beta*dxinv[0]*dxinv[0];
51 x(i,0,0,n) /= alpha*a(i,0,0) + dhx*(bX(i,0,0,n)+bX(i+1,0,0,n));
52}
53
54template <typename T>
56void mlabeclap_flux_x (Box const& box, Array4<T> const& fx, Array4<T const> const& sol,
57 Array4<T const> const& bx, T fac, int ncomp) noexcept
58{
59 const auto lo = amrex::lbound(box);
60 const auto hi = amrex::ubound(box);
61
62 for (int n = 0; n < ncomp; ++n) {
64 for (int i = lo.x; i <= hi.x; ++i) {
65 fx(i,0,0,n) = -fac*bx(i,0,0,n)*(sol(i,0,0,n)-sol(i-1,0,0,n));
66 }
67 }
68}
69
70template <typename T>
72void mlabeclap_flux_xface (Box const& box, Array4<T> const& fx, Array4<T const> const& sol,
73 Array4<T const> const& bx, T fac, int xlen, int ncomp) noexcept
74{
75 const auto lo = amrex::lbound(box);
76
77 for (int n = 0; n < ncomp; ++n) {
78 int i = lo.x;
79 fx(i,0,0,n) = -fac*bx(i,0,0,n)*(sol(i,0,0,n)-sol(i-1,0,0,n));
80 i += xlen;
81 fx(i,0,0,n) = -fac*bx(i,0,0,n)*(sol(i,0,0,n)-sol(i-1,0,0,n));
82 }
83}
84
85template <typename T>
87void abec_gsrb (int i, int, int, int n, Array4<T> const& phi, Array4<T const> const& rhs,
88 T alpha, Array4<T const> const& a,
89 T dhx,
90 Array4<T const> const& bX,
91 Array4<int const> const& m0,
92 Array4<int const> const& m1,
93 Array4<T const> const& f0,
94 Array4<T const> const& f1,
95 Box const& vbox, int redblack) noexcept
96{
97 if ((i+redblack)%2 == 0) {
98 const auto vlo = amrex::lbound(vbox);
99 const auto vhi = amrex::ubound(vbox);
100
101 T cf0 = (i == vlo.x && m0(vlo.x-1,0,0) > 0)
102 ? f0(vlo.x,0,0,n) : T(0.0);
103 T cf1 = (i == vhi.x && m1(vhi.x+1,0,0) > 0)
104 ? f1(vhi.x,0,0,n) : T(0.0);
105
106 T delta = dhx*(bX(i,0,0,n)*cf0 + bX(i+1,0,0,n)*cf1);
107
108 T gamma = alpha*a(i,0,0)
109 + dhx*( bX(i,0,0,n) + bX(i+1,0,0,n) );
110
111 T rho = dhx*(bX(i ,0 ,0,n)*phi(i-1,0 ,0,n)
112 + bX(i+1,0 ,0,n)*phi(i+1,0 ,0,n));
113
114 phi(i,0,0,n) = (rhs(i,0,0,n) + rho - phi(i,0,0,n)*delta)
115 / (gamma - delta);
116 }
117}
118
119template <typename T>
121void abec_gsrb_os (int i, int, int, int n, Array4<T> const& phi, Array4<T const> const& rhs,
122 T alpha, Array4<T const> const& a,
123 T dhx,
124 Array4<T const> const& bX,
125 Array4<int const> const& m0,
126 Array4<int const> const& m1,
127 Array4<T const> const& f0,
128 Array4<T const> const& f1,
129 Array4<int const> const& osm,
130 Box const& vbox, int redblack) noexcept
131{
132 if ((i+redblack)%2 == 0) {
133 if (osm(i,0,0) == 0) {
134 phi(i,0,0) = T(0.0);
135 } else {
136 const auto vlo = amrex::lbound(vbox);
137 const auto vhi = amrex::ubound(vbox);
138
139 T cf0 = (i == vlo.x && m0(vlo.x-1,0,0) > 0)
140 ? f0(vlo.x,0,0,n) : T(0.0);
141 T cf1 = (i == vhi.x && m1(vhi.x+1,0,0) > 0)
142 ? f1(vhi.x,0,0,n) : T(0.0);
143
144 T delta = dhx*(bX(i,0,0,n)*cf0 + bX(i+1,0,0,n)*cf1);
145
146 T gamma = alpha*a(i,0,0)
147 + dhx*( bX(i,0,0,n) + bX(i+1,0,0,n) );
148
149 T rho = dhx*(bX(i ,0 ,0,n)*phi(i-1,0 ,0,n)
150 + bX(i+1,0 ,0,n)*phi(i+1,0 ,0,n));
151
152 phi(i,0,0,n) = (rhs(i,0,0,n) + rho - phi(i,0,0,n)*delta)
153 / (gamma - delta);
154 }
155 }
156}
157
158template <typename T>
160void abec_jacobi (int i, int, int, int n, Array4<T> const& phi,
161 Array4<T const> const& rhs, Array4<T const> const& Ax,
162 T alpha, Array4<T const> const& a,
163 T dhx,
164 Array4<T const> const& bX,
165 Array4<int const> const& m0,
166 Array4<int const> const& m1,
167 Array4<T const> const& f0,
168 Array4<T const> const& f1,
169 Box const& vbox) noexcept
170{
171 const auto vlo = amrex::lbound(vbox);
172 const auto vhi = amrex::ubound(vbox);
173
174 T cf0 = (i == vlo.x && m0(vlo.x-1,0,0) > 0)
175 ? f0(vlo.x,0,0,n) : T(0.0);
176 T cf1 = (i == vhi.x && m1(vhi.x+1,0,0) > 0)
177 ? f1(vhi.x,0,0,n) : T(0.0);
178
179 T delta = dhx*(bX(i,0,0,n)*cf0 + bX(i+1,0,0,n)*cf1);
180
181 T gamma = alpha*a(i,0,0)
182 + dhx*( bX(i,0,0,n) + bX(i+1,0,0,n) );
183
184 phi(i,0,0,n) += T(2.0/3.0) * (rhs(i,0,0,n) - Ax(i,0,0,n)) / (gamma - delta);
185}
186
187template <typename T>
189void abec_jacobi_os (int i, int, int, int n, Array4<T> const& phi,
190 Array4<T const> const& rhs, Array4<T const> const& Ax,
191 T alpha, Array4<T const> const& a,
192 T dhx,
193 Array4<T const> const& bX,
194 Array4<int const> const& m0,
195 Array4<int const> const& m1,
196 Array4<T const> const& f0,
197 Array4<T const> const& f1,
198 Array4<int const> const& osm,
199 Box const& vbox) noexcept
200{
201 if (osm(i,0,0) == 0) {
202 phi(i,0,0) = T(0.0);
203 } else {
204 const auto vlo = amrex::lbound(vbox);
205 const auto vhi = amrex::ubound(vbox);
206
207 T cf0 = (i == vlo.x && m0(vlo.x-1,0,0) > 0)
208 ? f0(vlo.x,0,0,n) : T(0.0);
209 T cf1 = (i == vhi.x && m1(vhi.x+1,0,0) > 0)
210 ? f1(vhi.x,0,0,n) : T(0.0);
211
212 T delta = dhx*(bX(i,0,0,n)*cf0 + bX(i+1,0,0,n)*cf1);
213
214 T gamma = alpha*a(i,0,0)
215 + dhx*( bX(i,0,0,n) + bX(i+1,0,0,n) );
216
217 phi(i,0,0,n) += T(2.0/3.0) * (rhs(i,0,0,n) - Ax(i,0,0,n)) / (gamma - delta);
218 }
219}
220
221template <typename T>
224 Box const& /*box*/, Array4<T> const& /*phi*/, Array4<T const> const& /*rhs*/,
225 T /*alpha*/, Array4<T const> const& /*a*/,
226 T /*dhx*/,
227 Array4<T const> const& /*bX*/,
228 Array4<int const> const& /*m0*/,
229 Array4<int const> const& /*m1*/,
230 Array4<T const> const& /*f0*/,
231 Array4<T const> const& /*f1*/,
232 Box const& /*vbox*/, int /*redblack*/, int /*nc*/) noexcept
233{
234 amrex::Abort("abec_gsrb_with_line_solve not implemented in 1D");
235}
236
237template <typename T>
239void overset_rescale_bcoef_x (Box const& box, Array4<T> const& bX, Array4<int const> const& osm,
240 int ncomp, T osfac) noexcept
241{
242 const auto lo = amrex::lbound(box);
243 const auto hi = amrex::ubound(box);
244 for (int n = 0; n < ncomp; ++n) {
245 for (int i = lo.x; i <= hi.x; ++i) {
246 if ((osm(i-1,0,0)+osm(i,0,0)) == 1) {
247 bX(i,0,0,n) *= osfac;
248 }
249 }
250 }
251}
252
253}
254#endif
#define AMREX_PRAGMA_SIMD
Definition AMReX_Extension.H:80
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
Definition AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void overset_rescale_bcoef_x(Box const &box, Array4< T > const &bX, Array4< int const > const &osm, int ncomp, T osfac) noexcept
Definition AMReX_MLABecLap_1D_K.H:239
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_adotx_os(int i, int, int, int n, Array4< T > const &y, Array4< T const > const &x, Array4< T const > const &a, Array4< T const > const &bX, Array4< int const > const &osm, GpuArray< T, AMREX_SPACEDIM > const &dxinv, T alpha, T beta) noexcept
Definition AMReX_MLABecLap_1D_K.H:24
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_x(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, Array4< T const > const &bx, T fac, int ncomp) noexcept
Definition AMReX_MLABecLap_1D_K.H:56
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_gsrb_os(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Array4< int const > const &osm, Box const &vbox, int redblack) noexcept
Definition AMReX_MLABecLap_1D_K.H:121
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_jacobi(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Box const &vbox) noexcept
Definition AMReX_MLABecLap_1D_K.H:160
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_gsrb(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Box const &vbox, int redblack) noexcept
Definition AMReX_MLABecLap_1D_K.H:87
AMREX_FORCE_INLINE void abec_gsrb_with_line_solve(Box const &, Array4< T > const &, Array4< T const > const &, T, Array4< T const > const &, T, Array4< T const > const &, Array4< int const > const &, Array4< int const > const &, Array4< T const > const &, Array4< T const > const &, Box const &, int, int) noexcept
Definition AMReX_MLABecLap_1D_K.H:223
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_normalize(int i, int, int, int n, Array4< T > const &x, Array4< T const > const &a, Array4< T const > const &bX, GpuArray< T, AMREX_SPACEDIM > const &dxinv, T alpha, T beta) noexcept
Definition AMReX_MLABecLap_1D_K.H:44
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:225
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_jacobi_os(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Array4< int const > const &osm, Box const &vbox) noexcept
Definition AMReX_MLABecLap_1D_K.H:189
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_xface(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, Array4< T const > const &bx, T fac, int xlen, int ncomp) noexcept
Definition AMReX_MLABecLap_1D_K.H:72
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_adotx(int i, int, int, int n, Array4< T > const &y, Array4< T const > const &x, Array4< T const > const &a, Array4< T const > const &bX, GpuArray< T, AMREX_SPACEDIM > const &dxinv, T alpha, T beta) noexcept
Definition AMReX_MLABecLap_1D_K.H:9
Definition AMReX_Array4.H:61
Definition AMReX_Array.H:34