/home/runner/work/amr-wind/amr-wind/amr-wind/equation_systems/CompRHSOps.H Source File

AMR-Wind API: /home/runner/work/amr-wind/amr-wind/amr-wind/equation_systems/CompRHSOps.H Source File
AMR-Wind API v0.1.0
CFD solver for wind plant simulations
Loading...
Searching...
No Matches
CompRHSOps.H
Go to the documentation of this file.
1#ifndef COMPRHSOPS_H
2#define COMPRHSOPS_H
3
4#include "amr-wind/incflo_enums.H"
5#include "amr-wind/equation_systems/PDEOps.H"
6#include "amr-wind/equation_systems/SchemeTraits.H"
7#include "AMReX_REAL.H"
8
9using namespace amrex::literals;
10
11namespace amr_wind::pde {
12
19template <typename PDE, typename Scheme>
20struct ComputeRHSOp
21{
22 explicit ComputeRHSOp(PDEFields& fields_in)
23 : fields(fields_in), density(fields_in.repo.get_field("density"))
24 {}
25
32 void predictor_rhs(
33 const DiffusionType difftype, const amrex::Real dt, bool mesh_mapping)
34 {
35 amrex::Real factor = 0.0_rt;
36 switch (difftype) {
37 case DiffusionType::Explicit:
38 factor = 1.0_rt;
39 break;
40
41 case DiffusionType::Crank_Nicolson:
42 factor = 0.5_rt;
43 break;
44
45 case DiffusionType::Implicit:
46 factor = 0.0_rt;
47 break;
48
49 default:
50 amrex::Abort("Invalid diffusion type");
51 }
52
53 // Field states for diffusion and advection terms. In Godunov scheme
54 // these terms only have one state.
55 auto fstate = std::is_same<Scheme, fvm::Godunov>::value
56 ? FieldState::New
57 : FieldState::Old;
58
59 const int nlevels = fields.repo.num_active_levels();
60
61 // for RHS evaluation velocity field should be in stretched space
62 auto& field = fields.field;
63 if (field.in_uniform_space() && mesh_mapping) {
64 field.to_stretched_space();
65 }
66 auto& field_old = field.state(FieldState::Old);
67 if (field_old.in_uniform_space() && mesh_mapping) {
68 field_old.to_stretched_space();
69 }
70
71 auto& den_new = density.state(FieldState::New);
72 auto& den_old = density.state(FieldState::Old);
73 auto& src_term = fields.src_term;
74 auto& diff_term = fields.diff_term.state(fstate);
75 auto& conv_term = fields.conv_term.state(fstate);
76 auto& mask_cell = fields.repo.get_int_field("mask_cell");
77 Field const* mesh_detJ =
78 mesh_mapping ? &(fields.repo.get_mesh_mapping_det_j(FieldLoc::CELL))
79 : nullptr;
80
81 for (int lev = 0; lev < nlevels; ++lev) {
82 const auto& fld_arrs = field(lev).arrays();
83 const auto& fld_o_arrs = field_old(lev).const_arrays();
84 const auto& rho_o_arrs = den_old(lev).const_arrays();
85 const auto& rho_arrs = den_new(lev).const_arrays();
86 const auto& src_arrs = src_term(lev).const_arrays();
87 const auto& diff_arrs = diff_term(lev).const_arrays();
88 const auto& ddt_o_arrs = conv_term(lev).const_arrays();
89 const auto& imask_arrs = mask_cell(lev).const_arrays();
90 const auto& detJ_arrs =
91 mesh_mapping ? ((*mesh_detJ)(lev).const_arrays())
92 : amrex::MultiArray4<amrex::Real const>();
93
94 if (PDE::multiply_rho) {
95 // Remove multiplication by density as it will be added back
96 // in solver
97 amrex::ParallelFor(
98 field(lev), amrex::IntVect(0), PDE::ndim,
99 [=] AMREX_GPU_DEVICE(
100 int nbx, int i, int j, int k, int n) noexcept {
101 amrex::Real det_j =
102 mesh_mapping ? (detJ_arrs[nbx](i, j, k)) : 1.0_rt;
103
104 fld_arrs[nbx](i, j, k, n) =
105 rho_o_arrs[nbx](i, j, k) * det_j *
106 fld_o_arrs[nbx](i, j, k, n) +
107 static_cast<amrex::Real>(imask_arrs[nbx](i, j, k)) *
108 dt *
109 (ddt_o_arrs[nbx](i, j, k, n) +
110 det_j * src_arrs[nbx](i, j, k, n) +
111 factor * diff_arrs[nbx](i, j, k, n));
112
113 fld_arrs[nbx](i, j, k, n) /= rho_arrs[nbx](i, j, k);
114
115 if (difftype == DiffusionType::Explicit) {
116 fld_arrs[nbx](i, j, k, n) /= det_j;
117 }
118 });
119 } else {
120 amrex::ParallelFor(
121 field(lev), amrex::IntVect(0), PDE::ndim,
122 [=] AMREX_GPU_DEVICE(
123 int nbx, int i, int j, int k, int n) noexcept {
124 amrex::Real det_j =
125 mesh_mapping ? (detJ_arrs[nbx](i, j, k)) : 1.0_rt;
126
127 fld_arrs[nbx](i, j, k, n) =
128 det_j * fld_o_arrs[nbx](i, j, k, n) +
129 static_cast<amrex::Real>(imask_arrs[nbx](i, j, k)) *
130 dt *
131 (ddt_o_arrs[nbx](i, j, k, n) +
132 det_j * src_arrs[nbx](i, j, k, n) +
133 factor * diff_arrs[nbx](i, j, k, n));
134
135 if (difftype == DiffusionType::Explicit) {
136 fld_arrs[nbx](i, j, k, n) /= det_j;
137 }
138 });
139 }
140 }
141 amrex::Gpu::streamSynchronize();
142 }
143
150 void corrector_rhs(
151 const DiffusionType difftype, const amrex::Real dt, bool mesh_mapping)
152 {
153 amrex::Real ofac = 0.0_rt;
154 amrex::Real nfac = 0.0_rt;
155 switch (difftype) {
156 case DiffusionType::Explicit:
157 ofac = 0.5_rt;
158 nfac = 0.5_rt;
159 break;
160
161 case DiffusionType::Crank_Nicolson:
162 ofac = 0.5_rt;
163 nfac = 0.0_rt;
164 break;
165
166 case DiffusionType::Implicit:
167 ofac = 0.0_rt;
168 nfac = 0.0_rt;
169 break;
170
171 default:
172 amrex::Abort("Invalid diffusion type");
173 }
174
175 const int nlevels = fields.repo.num_active_levels();
176
177 // for RHS evaluation velocity field should be in stretched space
178 auto& field = fields.field;
179 if (field.in_uniform_space() && mesh_mapping) {
180 field.to_stretched_space();
181 }
182 auto& field_old = field.state(FieldState::Old);
183 if (field_old.in_uniform_space() && mesh_mapping) {
184 field_old.to_stretched_space();
185 }
186
187 auto& den_new = density.state(FieldState::New);
188 auto& den_old = density.state(FieldState::Old);
189 auto& src_term = fields.src_term;
190 auto& diff_term = fields.diff_term;
191 auto& conv_term = fields.conv_term;
192 auto& diff_term_old = fields.diff_term.state(FieldState::Old);
193 auto& conv_term_old = fields.conv_term.state(FieldState::Old);
194 auto& mask_cell = fields.repo.get_int_field("mask_cell");
195 Field const* mesh_detJ =
196 mesh_mapping ? &(fields.repo.get_mesh_mapping_det_j(FieldLoc::CELL))
197 : nullptr;
198
199 for (int lev = 0; lev < nlevels; ++lev) {
200 const auto& fld_arrs = field(lev).arrays();
201 const auto& fld_o_arrs = field_old(lev).const_arrays();
202 const auto& rho_o_arrs = den_old(lev).const_arrays();
203 const auto& rho_arrs = den_new(lev).const_arrays();
204 const auto& src_arrs = src_term(lev).const_arrays();
205 const auto& diff_arrs = diff_term(lev).const_arrays();
206 const auto& ddt_arrs = conv_term(lev).const_arrays();
207 const auto& diff_o_arrs = diff_term_old(lev).const_arrays();
208 const auto& ddt_o_arrs = conv_term_old(lev).const_arrays();
209 const auto& imask_arrs = mask_cell(lev).const_arrays();
210 const auto& detJ_arrs =
211 mesh_mapping ? ((*mesh_detJ)(lev).const_arrays())
212 : amrex::MultiArray4<amrex::Real const>();
213
214 if (PDE::multiply_rho) {
215 // Remove multiplication by density as it will be added back
216 // in solver
217 amrex::ParallelFor(
218 field(lev), amrex::IntVect(0), PDE::ndim,
219 [=] AMREX_GPU_DEVICE(
220 int nbx, int i, int j, int k, int n) noexcept {
221 amrex::Real det_j =
222 mesh_mapping ? (detJ_arrs[nbx](i, j, k)) : 1.0_rt;
223
224 fld_arrs[nbx](i, j, k, n) =
225 rho_o_arrs[nbx](i, j, k) * det_j *
226 fld_o_arrs[nbx](i, j, k, n) +
227 static_cast<amrex::Real>(imask_arrs[nbx](i, j, k)) *
228 dt *
229 (0.5_rt * (ddt_o_arrs[nbx](i, j, k, n) +
230 ddt_arrs[nbx](i, j, k, n)) +
231 ofac * diff_o_arrs[nbx](i, j, k, n) +
232 nfac * diff_arrs[nbx](i, j, k, n) +
233 det_j * src_arrs[nbx](i, j, k, n));
234
235 fld_arrs[nbx](i, j, k, n) /= rho_arrs[nbx](i, j, k);
236
237 if (difftype == DiffusionType::Explicit) {
238 fld_arrs[nbx](i, j, k, n) /= det_j;
239 }
240 });
241 } else {
242 amrex::ParallelFor(
243 field(lev), amrex::IntVect(0), PDE::ndim,
244 [=] AMREX_GPU_DEVICE(
245 int nbx, int i, int j, int k, int n) noexcept {
246 amrex::Real det_j =
247 mesh_mapping ? (detJ_arrs[nbx](i, j, k)) : 1.0_rt;
248
249 fld_arrs[nbx](i, j, k, n) =
250 det_j * fld_o_arrs[nbx](i, j, k, n) +
251 static_cast<amrex::Real>(imask_arrs[nbx](i, j, k)) *
252 dt *
253 (0.5_rt * (ddt_o_arrs[nbx](i, j, k, n) +
254 ddt_arrs[nbx](i, j, k, n)) +
255 ofac * diff_o_arrs[nbx](i, j, k, n) +
256 nfac * diff_arrs[nbx](i, j, k, n) +
257 det_j * src_arrs[nbx](i, j, k, n));
258
259 if (difftype == DiffusionType::Explicit) {
260 fld_arrs[nbx](i, j, k, n) /= det_j;
261 }
262 });
263 }
264 }
265 amrex::Gpu::streamSynchronize();
266 }
267
268 void improve_explicit_diff(const amrex::Real dt)
269 {
270 explicit_rk2_post_hoc(dt);
271 }
272
273 void explicit_rk2_post_hoc(const amrex::Real dt)
274 {
275 // divtau must be the difference between new and old before calling this
276 const auto& d_divtau = fields.diff_term;
277
278 auto& dof = fields.field;
279 const auto& repo = fields.repo;
280 const auto& mask_cell = repo.get_int_field("mask_cell");
281
282 const int nlevels = repo.num_active_levels();
283 for (int lev = 0; lev < nlevels; ++lev) {
284 auto f_arrs = dof(lev).arrays();
285 const auto& d_diff_arrs = d_divtau(lev).const_arrays();
286 const auto& mask_arrs = mask_cell(lev).const_arrays();
287 const auto& rho_arrs = density(lev).const_arrays();
288
289 amrex::ParallelFor(
290 dof(lev), amrex::IntVect(0), dof.num_comp(),
291 [=] AMREX_GPU_DEVICE(
292 int nbx, int i, int j, int k, int n) noexcept {
293 auto factor =
294 0.5_rt * dt * (amrex::Real)mask_arrs[nbx](i, j, k);
295 if (PDE::multiply_rho) {
296 factor /= rho_arrs[nbx](i, j, k);
297 }
298 f_arrs[nbx](i, j, k, n) +=
299 factor * d_diff_arrs[nbx](i, j, k, n);
300 });
301 }
302 }
303
304 // data members
305 PDEFields& fields;
306 Field& density;
307};
308
309} // namespace amr_wind::pde
310
311#endif /* COMPRHSOPS_H */
amr_wind::Field
Definition Field.H:116
amr_wind::FieldLoc::CELL
@ CELL
Cell-centered (default)
Definition FieldDescTypes.H:28
amr_wind::FieldState::New
@ New
Same as FieldState::NP1.
Definition FieldDescTypes.H:20
amr_wind::FieldState::Old
@ Old
Same as FieldState::N.
Definition FieldDescTypes.H:21
incflo_enums.H
DiffusionType
DiffusionType
Definition incflo_enums.H:4
DiffusionType::Implicit
@ Implicit
Definition incflo_enums.H:4
DiffusionType::Crank_Nicolson
@ Crank_Nicolson
Definition incflo_enums.H:4
DiffusionType::Explicit
@ Explicit
Definition incflo_enums.H:4
amr_wind::pde
Definition AdvOp_Godunov.H:21
amr_wind::pde::ComputeRHSOp::ComputeRHSOp
ComputeRHSOp(PDEFields &fields_in)
Definition CompRHSOps.H:22
amr_wind::pde::ComputeRHSOp::predictor_rhs
void predictor_rhs(const DiffusionType difftype, const amrex::Real dt, bool mesh_mapping)
Definition CompRHSOps.H:32
amr_wind::pde::ComputeRHSOp::improve_explicit_diff
void improve_explicit_diff(const amrex::Real dt)
Definition CompRHSOps.H:268
amr_wind::pde::ComputeRHSOp::explicit_rk2_post_hoc
void explicit_rk2_post_hoc(const amrex::Real dt)
Definition CompRHSOps.H:273
amr_wind::pde::ComputeRHSOp::corrector_rhs
void corrector_rhs(const DiffusionType difftype, const amrex::Real dt, bool mesh_mapping)
Definition CompRHSOps.H:150
amr_wind::pde::ComputeRHSOp::fields
PDEFields & fields
Definition CompRHSOps.H:305
amr_wind::pde::ComputeRHSOp::density
Field & density
Definition CompRHSOps.H:306
amr_wind::pde::PDEFields
Definition PDEFields.H:27
Generated by doxygen 1.14.0