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