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