/home/runner/work/amr-wind/amr-wind/amr-wind/overset/overset_ops_K.H Source File

AMR-Wind API: /home/runner/work/amr-wind/amr-wind/amr-wind/overset/overset_ops_K.H Source File
AMR-Wind API v0.1.0
CFD solver for wind plant simulations
Loading...
Searching...
No Matches
overset_ops_K.H
Go to the documentation of this file.
1#ifndef OVERSET_OPS_K_H_
2#define OVERSET_OPS_K_H_
3
4#include "amr-wind/core/vs/vector_space.H"
5#include <AMReX_FArrayBox.H>
6#include "AMReX_REAL.H"
7
8using namespace amrex::literals;
9
10namespace amr_wind::overset_ops {
11
12// Approximate signed distance function
13amrex::Real AMREX_GPU_DEVICE AMREX_FORCE_INLINE
14asdf(const amrex::Real a_vof, const amrex::Real i_th, const amrex::Real tiny)
15{
16 // function of local vof value and interface thickness
17 return (i_th * std::log((a_vof + tiny) / (1.0_rt - a_vof + tiny)));
18}
19
20amrex::Real AMREX_GPU_DEVICE AMREX_FORCE_INLINE alpha_flux(
21 const int i,
22 const int j,
23 const int k,
24 const int dir,
25 const amrex::Real margin,
26 amrex::Array4<amrex::Real const> const& vof,
27 amrex::Array4<amrex::Real const> const& tg_vof,
28 amrex::Array4<amrex::Real const> const& normal)
29{
30 // Set up neighbor indices
31 const amrex::IntVect iv{i, j, k};
32 const amrex::IntVect dv{(int)(dir == 0), (int)(dir == 1), (int)(dir == 2)};
33 const amrex::IntVect ivm = iv - dv;
34
35 // Gradient of phi normal to interface
36 const amrex::Real gphi = (vof(iv) - vof(ivm));
37 // Normal vector in each cell (already normalized)
38 const amrex::Real norm_ = std::abs(normal(iv, dir));
39 const amrex::Real norm_nb = std::abs(normal(ivm, dir));
40
41 // Determine which delta_phi (and multiply by normal)
42 // The sign depends on side of flux face
43 const amrex::Real dphi_ = (vof(iv) - tg_vof(iv)) * norm_;
44 const amrex::Real dphi_nb = (tg_vof(ivm) - vof(ivm)) * norm_nb;
45 // Average value used across the interface
46 amrex::Real dphi_eval = 0.5_rt * (dphi_ + dphi_nb);
47 // Upwinding when on the gas side, downwinding on the liquid
48 // Across the interface defined as crossing 0.5_rt or within margin of
49 // 0.5_rt
50 if ((std::abs(vof(iv) - 0.5_rt) > margin ||
51 std::abs(vof(ivm) - 0.5_rt) > margin)) {
52 if (gphi > 0.0_rt) {
53 dphi_eval = (vof(ivm) < 0.5_rt && vof(iv) <= 0.5_rt + margin)
54 ? dphi_nb
55 : dphi_eval;
56 dphi_eval = (vof(ivm) >= 0.5_rt - margin && vof(iv) > 0.5_rt)
57 ? dphi_
58 : dphi_eval;
59 }
60 if (gphi < 0.0_rt) {
61 dphi_eval = (vof(iv) < 0.5_rt && vof(ivm) <= 0.5_rt + margin)
62 ? dphi_
63 : dphi_eval;
64 dphi_eval = (vof(iv) >= 0.5_rt - margin && vof(ivm) > 0.5_rt)
65 ? dphi_nb
66 : dphi_eval;
67 }
68 }
69 return dphi_eval;
70}
71
72void AMREX_GPU_DEVICE AMREX_FORCE_INLINE velocity_face(
73 const int i,
74 const int j,
75 const int k,
76 const int dir,
77 amrex::Array4<amrex::Real const> const& vof,
78 amrex::Array4<amrex::Real const> const& velocity,
79 amrex::Real& uface,
80 amrex::Real& vface,
81 amrex::Real& wface)
82{
83 // Set up neighbor indices
84 const amrex::IntVect iv{i, j, k};
85 const amrex::IntVect dv{(int)(dir == 0), (int)(dir == 1), (int)(dir == 2)};
86 const amrex::IntVect ivm = iv - dv;
87
88 // Gradient of phi normal to interface
89 const amrex::Real gphi = (vof(iv) - vof(ivm));
90
91 // Get velocities on both sides
92 const amrex::Real u_ = velocity(iv, 0);
93 const amrex::Real v_ = velocity(iv, 1);
94 const amrex::Real w_ = velocity(iv, 2);
95 const amrex::Real u_nb = velocity(ivm, 0);
96 const amrex::Real v_nb = velocity(ivm, 1);
97 const amrex::Real w_nb = velocity(ivm, 2);
98 // Average value when gphi = 0
99 uface = 0.5_rt * (u_ + u_nb);
100 vface = 0.5_rt * (v_ + v_nb);
101 wface = 0.5_rt * (w_ + w_nb);
102 // Use simple upwinding
103 if (gphi > 0.0_rt) {
104 uface = u_nb;
105 vface = v_nb;
106 wface = w_nb;
107 }
108 if (gphi < 0.0_rt) {
109 uface = u_;
110 vface = v_;
111 wface = w_;
112 }
113}
114
115void AMREX_GPU_DEVICE AMREX_FORCE_INLINE gp_rho_face(
116 const int i,
117 const int j,
118 const int k,
119 const int dir,
120 amrex::Array4<amrex::Real const> const& vof,
121 amrex::Array4<amrex::Real const> const& gp,
122 amrex::Array4<amrex::Real const> const& rho,
123 amrex::Real& uface,
124 amrex::Real& vface,
125 amrex::Real& wface)
126{
127 // Set up neighbor indices
128 const amrex::IntVect iv{i, j, k};
129 const amrex::IntVect dv{(int)(dir == 0), (int)(dir == 1), (int)(dir == 2)};
130 const amrex::IntVect ivm = iv - dv;
131
132 // Gradient of phi normal to interface
133 const amrex::Real gphi = (vof(iv) - vof(ivm));
134
135 // Get velocities on both sides
136 const amrex::Real u_ = gp(iv, 0) / rho(iv);
137 const amrex::Real v_ = gp(iv, 1) / rho(iv);
138 const amrex::Real w_ = gp(iv, 2) / rho(iv);
139 const amrex::Real u_nb = gp(ivm, 0) / rho(ivm);
140 const amrex::Real v_nb = gp(ivm, 1) / rho(ivm);
141 const amrex::Real w_nb = gp(ivm, 2) / rho(ivm);
142 // Average value when gphi = 0
143 uface = 0.5_rt * (u_ + u_nb);
144 vface = 0.5_rt * (v_ + v_nb);
145 wface = 0.5_rt * (w_ + w_nb);
146 // Use simple upwinding
147 if (gphi > 0.0_rt) {
148 uface = u_nb;
149 vface = v_nb;
150 wface = w_nb;
151 }
152 if (gphi < 0.0_rt) {
153 uface = u_;
154 vface = v_;
155 wface = w_;
156 }
157}
158
159vs::Tensor AMREX_GPU_DEVICE AMREX_FORCE_INLINE gp_flux_tensor(
160 const int i,
161 const int j,
162 const int k,
163 amrex::Array4<amrex::Real const> const& fx,
164 amrex::Array4<amrex::Real const> const& fy,
165 amrex::Array4<amrex::Real const> const& fz,
166 const amrex::Real tiny)
167{
168 // Averaging depends on iblank array, encoded in 8th component
169 const amrex::Real avg_fx = fx(i, j, k, 8) + fx(i, j, k - 1, 8) +
170 fx(i, j - 1, k, 8) + fx(i, j - 1, k - 1, 8) +
171 tiny;
172 const amrex::Real avg_fy = fy(i, j, k, 8) + fy(i - 1, j, k, 8) +
173 fy(i, j, k - 1, 8) + fy(i - 1, j, k - 1, 8) +
174 tiny;
175 const amrex::Real avg_fz = fz(i, j, k, 8) + fz(i - 1, j, k, 8) +
176 fz(i, j - 1, k, 8) + fz(i - 1, j - 1, k, 8) +
177 tiny;
178
179 auto f_otimes_gradp = vs::Tensor::identity();
180 // Average fluxes (faces) to nodes where pressure exists
181 f_otimes_gradp.xx() = (fx(i, j, k, 5) + fx(i, j, k - 1, 5) +
182 fx(i, j - 1, k, 5) + fx(i, j - 1, k - 1, 5)) /
183 avg_fx;
184 f_otimes_gradp.xy() = (fx(i, j, k, 6) + fx(i, j, k - 1, 6) +
185 fx(i, j - 1, k, 6) + fx(i, j - 1, k - 1, 6)) /
186 avg_fx;
187 f_otimes_gradp.xz() = (fx(i, j, k, 7) + fx(i, j, k - 1, 7) +
188 fx(i, j - 1, k, 7) + fx(i, j - 1, k - 1, 7)) /
189 avg_fx;
190 f_otimes_gradp.yx() = (fy(i, j, k, 5) + fy(i - 1, j, k, 5) +
191 fy(i, j, k - 1, 5) + fy(i - 1, j, k - 1, 5)) /
192 avg_fy;
193 f_otimes_gradp.yy() = (fy(i, j, k, 6) + fy(i - 1, j, k, 6) +
194 fy(i, j, k - 1, 6) + fy(i - 1, j, k - 1, 6)) /
195 avg_fy;
196 f_otimes_gradp.yz() = (fy(i, j, k, 7) + fy(i - 1, j, k, 7) +
197 fy(i, j, k - 1, 7) + fy(i - 1, j, k - 1, 7)) /
198 avg_fy;
199 f_otimes_gradp.zx() = (fz(i, j, k, 5) + fz(i - 1, j, k, 5) +
200 fz(i, j - 1, k, 5) + fz(i - 1, j - 1, k, 5)) /
201 avg_fz;
202 f_otimes_gradp.zy() = (fz(i, j, k, 6) + fz(i - 1, j, k, 6) +
203 fz(i, j - 1, k, 6) + fz(i - 1, j - 1, k, 6)) /
204 avg_fz;
205 f_otimes_gradp.zz() = (fz(i, j, k, 7) + fz(i - 1, j, k, 7) +
206 fz(i, j - 1, k, 7) + fz(i - 1, j - 1, k, 7)) /
207 avg_fz;
208
209 return f_otimes_gradp;
210}
211
212vs::Tensor AMREX_GPU_DEVICE AMREX_FORCE_INLINE normal_reinit_tensor(
213 const int i,
214 const int j,
215 const int k,
216 amrex::Array4<amrex::Real const> const& fx,
217 amrex::Array4<amrex::Real const> const& fy,
218 amrex::Array4<amrex::Real const> const& fz,
219 amrex::Array4<amrex::Real const> const& vof,
220 const amrex::Real tiny)
221{
222 // Averaging depends on iblank array, encoded in 8th component
223 const amrex::Real avg_fx = fx(i, j, k, 8) + fx(i, j, k - 1, 8) +
224 fx(i, j - 1, k, 8) + fx(i, j - 1, k - 1, 8) +
225 tiny;
226 const amrex::Real avg_fy = fy(i, j, k, 8) + fy(i - 1, j, k, 8) +
227 fy(i, j, k - 1, 8) + fy(i - 1, j, k - 1, 8) +
228 tiny;
229 const amrex::Real avg_fz = fz(i, j, k, 8) + fz(i - 1, j, k, 8) +
230 fz(i, j - 1, k, 8) + fz(i - 1, j - 1, k, 8) +
231 tiny;
232
233 auto n_zeta = vs::Vector::one();
234 // Get components of normal in each position
235 n_zeta.x() =
236 ((vof(i, j, k) - vof(i - 1, j, k)) * fx(i, j, k, 8) +
237 (vof(i, j - 1, k) - vof(i - 1, j - 1, k)) * fx(i, j - 1, k, 8) +
238 (vof(i, j, k - 1) - vof(i - 1, j, k - 1)) * fx(i, j, k - 1, 8) +
239 (vof(i, j - 1, k - 1) - vof(i - 1, j - 1, k - 1)) *
240 fx(i, j - 1, k - 1, 8)) /
241 avg_fx;
242 n_zeta.y() =
243 ((vof(i, j, k) - vof(i, j - 1, k)) * fy(i, j, k, 8) +
244 (vof(i - 1, j, k) - vof(i - 1, j - 1, k)) * fy(i - 1, j, k, 8) +
245 (vof(i, j, k - 1) - vof(i, j - 1, k - 1)) * fy(i, j, k - 1, 8) +
246 (vof(i - 1, j, k - 1) - vof(i - 1, j - 1, k - 1)) *
247 fy(i - 1, j, k - 1, 8)) /
248 avg_fy;
249 n_zeta.z() =
250 ((vof(i, j, k) - vof(i, j, k - 1)) * fz(i, j, k, 8) +
251 (vof(i - 1, j, k) - vof(i - 1, j, k - 1)) * fz(i - 1, j, k, 8) +
252 (vof(i, j - 1, k) - vof(i, j - 1, k - 1)) * fz(i, j - 1, k, 8) +
253 (vof(i - 1, j - 1, k) - vof(i - 1, j - 1, k - 1)) *
254 fz(i - 1, j - 1, k, 8)) /
255 avg_fz;
256
257 n_zeta.normalize();
258
259 // To do outer product between vectors, set up tensors
260 auto n_tensor = vs::Tensor::identity();
261 auto n_tensor_T = vs::Tensor::identity();
262 n_tensor.rows(n_zeta, n_zeta, n_zeta);
263 n_tensor_T.cols(n_zeta, n_zeta, n_zeta);
264 // Multiply tensors elementwise
265 auto n_otimes_n = vs::Tensor::identity();
266 n_otimes_n.rows(
267 n_tensor.x() * n_tensor_T.x(), n_tensor.y() * n_tensor_T.y(),
268 n_tensor.z() * n_tensor_T.z());
269
270 return n_otimes_n;
271}
272
273} // namespace amr_wind::overset_ops
274
275#endif
amr_wind::overset_ops
Definition overset_ops_K.H:10
amr_wind::overset_ops::alpha_flux
amrex::Real AMREX_GPU_DEVICE AMREX_FORCE_INLINE alpha_flux(const int i, const int j, const int k, const int dir, const amrex::Real margin, amrex::Array4< amrex::Real const > const &vof, amrex::Array4< amrex::Real const > const &tg_vof, amrex::Array4< amrex::Real const > const &normal)
Definition overset_ops_K.H:20
amr_wind::overset_ops::velocity_face
void AMREX_GPU_DEVICE AMREX_FORCE_INLINE velocity_face(const int i, const int j, const int k, const int dir, amrex::Array4< amrex::Real const > const &vof, amrex::Array4< amrex::Real const > const &velocity, amrex::Real &uface, amrex::Real &vface, amrex::Real &wface)
Definition overset_ops_K.H:72
amr_wind::overset_ops::gp_flux_tensor
vs::Tensor AMREX_GPU_DEVICE AMREX_FORCE_INLINE gp_flux_tensor(const int i, const int j, const int k, amrex::Array4< amrex::Real const > const &fx, amrex::Array4< amrex::Real const > const &fy, amrex::Array4< amrex::Real const > const &fz, const amrex::Real tiny)
Definition overset_ops_K.H:159
amr_wind::overset_ops::asdf
amrex::Real AMREX_GPU_DEVICE AMREX_FORCE_INLINE asdf(const amrex::Real a_vof, const amrex::Real i_th, const amrex::Real tiny)
Definition overset_ops_K.H:14
amr_wind::overset_ops::normal_reinit_tensor
vs::Tensor AMREX_GPU_DEVICE AMREX_FORCE_INLINE normal_reinit_tensor(const int i, const int j, const int k, amrex::Array4< amrex::Real const > const &fx, amrex::Array4< amrex::Real const > const &fy, amrex::Array4< amrex::Real const > const &fz, amrex::Array4< amrex::Real const > const &vof, const amrex::Real tiny)
Definition overset_ops_K.H:212
amr_wind::overset_ops::gp_rho_face
void AMREX_GPU_DEVICE AMREX_FORCE_INLINE gp_rho_face(const int i, const int j, const int k, const int dir, amrex::Array4< amrex::Real const > const &vof, amrex::Array4< amrex::Real const > const &gp, amrex::Array4< amrex::Real const > const &rho, amrex::Real &uface, amrex::Real &vface, amrex::Real &wface)
Definition overset_ops_K.H:115
amr_wind::vs::Tensor
TensorT< amrex::Real > Tensor
Definition tensor.H:189
amr_wind::vs::TensorT< amrex::Real >::identity
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr TensorT< amrex::Real > identity() noexcept
Definition tensor.H:65
amr_wind::vs::VectorT< amrex::Real >::one
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr VectorT< amrex::Real > one()
Definition vector.H:50
vector_space.H
Generated by doxygen 1.14.0