/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#ifdef AMREX_USE_OMP
80#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
81#endif
82 for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
83 const auto& bx = mfi.tilebox();
84 auto fld = field(lev).array(mfi);
85 const auto fld_o = field_old(lev).const_array(mfi);
86 const auto rho_o = den_old(lev).const_array(mfi);
87 const auto rho = den_new(lev).const_array(mfi);
88 const auto src = src_term(lev).const_array(mfi);
89 const auto diff = diff_term(lev).const_array(mfi);
90 const auto ddt_o = conv_term(lev).const_array(mfi);
91 const auto imask = mask_cell(lev).const_array(mfi);
92 amrex::Array4<amrex::Real const> detJ =
93 mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi))
94 : amrex::Array4<amrex::Real const>();
95
96 if (PDE::multiply_rho) {
97 // Remove multiplication by density as it will be added back
98 // in solver
99 amrex::ParallelFor(
100 bx, PDE::ndim,
101 [=] AMREX_GPU_DEVICE(
102 int i, int j, int k, int n) noexcept {
103 amrex::Real det_j =
104 mesh_mapping ? (detJ(i, j, k)) : 1.0;
105
106 fld(i, j, k, n) =
107 rho_o(i, j, k) * det_j * fld_o(i, j, k, n) +
108 static_cast<amrex::Real>(imask(i, j, k)) * dt *
109 (ddt_o(i, j, k, n) +
110 det_j * src(i, j, k, n) +
111 factor * diff(i, j, k, n));
112
113 fld(i, j, k, n) /= rho(i, j, k);
114
115 if (difftype == DiffusionType::Explicit) {
116 fld(i, j, k, n) /= det_j;
117 }
118 });
119 } else {
120 amrex::ParallelFor(
121 bx, PDE::ndim,
122 [=] AMREX_GPU_DEVICE(
123 int i, int j, int k, int n) noexcept {
124 amrex::Real det_j =
125 mesh_mapping ? (detJ(i, j, k)) : 1.0;
126
127 fld(i, j, k, n) =
128 det_j * fld_o(i, j, k, n) +
129 static_cast<amrex::Real>(imask(i, j, k)) * dt *
130 (ddt_o(i, j, k, n) +
131 det_j * src(i, j, k, n) +
132 factor * diff(i, j, k, n));
133
134 if (difftype == DiffusionType::Explicit) {
135 fld(i, j, k, n) /= det_j;
136 }
137 });
138 }
139 }
140 }
141 }
142
149 void corrector_rhs(
150 const DiffusionType difftype, const amrex::Real dt, bool mesh_mapping)
151 {
152 amrex::Real ofac = 0.0;
153 amrex::Real nfac = 0.0;
154 switch (difftype) {
155 case DiffusionType::Explicit:
156 ofac = 0.5;
157 nfac = 0.5;
158 break;
159
160 case DiffusionType::Crank_Nicolson:
161 ofac = 0.5;
162 nfac = 0.0;
163 break;
164
165 case DiffusionType::Implicit:
166 ofac = 0.0;
167 nfac = 0.0;
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#ifdef AMREX_USE_OMP
200#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
201#endif
202 for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
203 const auto& bx = mfi.tilebox();
204 auto fld = field(lev).array(mfi);
205 const auto fld_o = field_old(lev).const_array(mfi);
206 const auto rho_o = den_old(lev).const_array(mfi);
207 const auto rho = den_new(lev).const_array(mfi);
208 const auto src = src_term(lev).const_array(mfi);
209 const auto diff = diff_term(lev).const_array(mfi);
210 const auto ddt = conv_term(lev).const_array(mfi);
211 const auto diff_o = diff_term_old(lev).const_array(mfi);
212 const auto ddt_o = conv_term_old(lev).const_array(mfi);
213 const auto imask = mask_cell(lev).const_array(mfi);
214 amrex::Array4<amrex::Real const> detJ =
215 mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi))
216 : amrex::Array4<amrex::Real const>();
217
218 if (PDE::multiply_rho) {
219 // Remove multiplication by density as it will be added back
220 // in solver
221 amrex::ParallelFor(
222 bx, PDE::ndim,
223 [=] AMREX_GPU_DEVICE(
224 int i, int j, int k, int n) noexcept {
225 amrex::Real det_j =
226 mesh_mapping ? (detJ(i, j, k)) : 1.0;
227
228 fld(i, j, k, n) =
229 rho_o(i, j, k) * det_j * fld_o(i, j, k, n) +
230 static_cast<amrex::Real>(imask(i, j, k)) * dt *
231 (0.5 *
232 (ddt_o(i, j, k, n) + ddt(i, j, k, n)) +
233 ofac * diff_o(i, j, k, n) +
234 nfac * diff(i, j, k, n) +
235 det_j * src(i, j, k, n));
236
237 fld(i, j, k, n) /= rho(i, j, k);
238
239 if (difftype == DiffusionType::Explicit) {
240 fld(i, j, k, n) /= det_j;
241 }
242 });
243 } else {
244 amrex::ParallelFor(
245 bx, PDE::ndim,
246 [=] AMREX_GPU_DEVICE(
247 int i, int j, int k, int n) noexcept {
248 amrex::Real det_j =
249 mesh_mapping ? (detJ(i, j, k)) : 1.0;
250
251 fld(i, j, k, n) =
252 det_j * fld_o(i, j, k, n) +
253 static_cast<amrex::Real>(imask(i, j, k)) * dt *
254 (0.5 *
255 (ddt_o(i, j, k, n) + ddt(i, j, k, n)) +
256 ofac * diff_o(i, j, k, n) +
257 nfac * diff(i, j, k, n) +
258 det_j * src(i, j, k, n));
259
260 if (difftype == DiffusionType::Explicit) {
261 fld(i, j, k, n) /= det_j;
262 }
263 });
264 }
265 }
266 }
267 }
268
269 // data members
270 PDEFields& fields;
271 Field& density;
272};
273
274} // namespace amr_wind::pde
275
276#endif /* COMPRHSOPS_H */
amr_wind::Field
Definition Field.H:116
amr_wind::Field::state
Field & state(const FieldState fstate)
Return field at a different time state.
Definition Field.cpp:114
amr_wind::FieldRepo::num_active_levels
int num_active_levels() const noexcept
Total number of levels currently active in the AMR mesh.
Definition FieldRepo.H:361
amr_wind::FieldRepo::get_mesh_mapping_det_j
Field & get_mesh_mapping_det_j(FieldLoc floc) const
Definition FieldRepo.cpp:190
amr_wind::FieldRepo::get_int_field
IntField & get_int_field(const std::string &name, const FieldState fstate=FieldState::New) const
Return a reference to an integer field.
Definition FieldRepo.cpp:276
amr_wind::FieldLoc::CELL
@ CELL
Cell-centered (default)
amr_wind::FieldState::New
@ New
Same as FieldState::NP1.
amr_wind::FieldState::Old
@ Old
Same as FieldState::N.
incflo_enums.H
DiffusionType
DiffusionType
Definition incflo_enums.H:4
DiffusionType::Implicit
@ Implicit
DiffusionType::Crank_Nicolson
@ Crank_Nicolson
DiffusionType::Explicit
@ Explicit
amr_wind::pde
Definition AdvOp_Godunov.H:16
amr_wind::pde::ComputeRHSOp
Definition CompRHSOps.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::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:270
amr_wind::pde::ComputeRHSOp::density
Field & density
Definition CompRHSOps.H:271
amr_wind::pde::PDEFields
Definition PDEFields.H:27
amr_wind::pde::PDEFields::src_term
Field & src_term
Source term for this PDE.
Definition PDEFields.H:39
amr_wind::pde::PDEFields::diff_term
Field & diff_term
Diffusion term for this PDE.
Definition PDEFields.H:41
amr_wind::pde::PDEFields::conv_term
Field & conv_term
Convective term for this PDE.
Definition PDEFields.H:43
amr_wind::pde::PDEFields::repo
FieldRepo & repo
Reference to the field repository instance.
Definition PDEFields.H:31
amr_wind::pde::PDEFields::field
Field & field
Solution variable (e.g., velocity, temperature)
Definition PDEFields.H:34
Generated by doxygen 1.12.0