/home/runner/work/kynema-sgf/kynema-sgf/src/fvm/curvature.H Source File

Kynema-SGF API: /home/runner/work/kynema-sgf/kynema-sgf/src/fvm/curvature.H Source File
Kynema-SGF API v0.1.0
CFD solver for wind plant simulations
Loading...
Searching...
No Matches
curvature.H
Go to the documentation of this file.
1#ifndef CURVATURE_H
2#define CURVATURE_H
3
4#include "src/fvm/fvm_utils.H"
5#include "AMReX_REAL.H"
6
7using namespace amrex::literals;
8
9namespace kynema_sgf::fvm {
10
14template <typename FTypeIn, typename FTypeOut>
15struct Curvature
16{
17 Curvature(FTypeOut& curphi, const FTypeIn& phi)
18 : m_curphi(curphi), m_phi(phi)
19 {
20 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
21 m_phi.num_comp() == m_curphi.num_comp(),
22 "m_phi.num_comp() = " + std::to_string(m_phi.num_comp()) +
23 ", m_curphi.num_comp() = " +
24 std::to_string(m_curphi.num_comp()));
25 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
26 m_phi.num_grow() > amrex::IntVect(0),
27 "m_phi.num_grow() = (" + std::to_string(m_phi.num_grow()[0]) +
28 ", " + std::to_string(m_phi.num_grow()[1]) + ", " +
29 std::to_string(m_phi.num_grow()[2]) + ")");
30 }
31
32 template <typename Stencil>
33 void apply(const int lev) const
34 {
35 const auto& geom = m_phi.repo().mesh().Geom(lev);
36 const auto& idx = geom.InvCellSizeArray();
37 const auto dlo = geom.Domain().smallEnd();
38 const auto dhi = geom.Domain().bigEnd();
39 const int ncomp = m_phi.num_comp();
40
41 auto& out_mf = m_curphi(lev);
42 auto const& in_arrs = m_phi(lev).const_arrays();
43 auto const& out_arrs = out_mf.arrays();
44
45 amrex::ParallelFor(
46 out_mf, amrex::IntVect(0),
47 [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) {
48 if (!Stencil::applies_to(i, j, k, dlo, dhi)) {
49 return;
50 }
51 auto const& phi = in_arrs[nbx];
52 auto const& curphi = out_arrs[nbx];
53 for (int n = 0; n < ncomp; ++n) {
54 const amrex::Real phixx =
55 ((Stencil::s00 * phi(i + 1, j, k, n)) +
56 (Stencil::s01 * phi(i, j, k, n)) +
57 (Stencil::s02 * phi(i - 1, j, k, n))) *
58 idx[0] * idx[0];
59 const amrex::Real phiyy =
60 ((Stencil::s10 * phi(i, j + 1, k, n)) +
61 (Stencil::s11 * phi(i, j, k, n)) +
62 (Stencil::s12 * phi(i, j - 1, k, n))) *
63 idx[1] * idx[1];
64 const amrex::Real phizz =
65 ((Stencil::s20 * phi(i, j, k + 1, n)) +
66 (Stencil::s21 * phi(i, j, k, n)) +
67 (Stencil::s22 * phi(i, j, k - 1, n))) *
68 idx[2] * idx[2];
69 const amrex::Real phiz =
70 ((Stencil::c20 * phi(i, j, k + 1, n)) +
71 (Stencil::c21 * phi(i, j, k, n)) +
72 (Stencil::c22 * phi(i, j, k - 1, n))) *
73 idx[2];
74 const amrex::Real phiz_ip1 =
75 ((Stencil::c20 * phi(i + 1, j, k + 1, n)) +
76 (Stencil::c21 * phi(i + 1, j, k, n)) +
77 (Stencil::c22 * phi(i + 1, j, k - 1, n))) *
78 idx[2];
79 const amrex::Real phiz_im1 =
80 ((Stencil::c20 * phi(i - 1, j, k + 1, n)) +
81 (Stencil::c21 * phi(i - 1, j, k, n)) +
82 (Stencil::c22 * phi(i - 1, j, k - 1, n))) *
83 idx[2];
84 const amrex::Real phiz_jp1 =
85 ((Stencil::c20 * phi(i, j + 1, k + 1, n)) +
86 (Stencil::c21 * phi(i, j + 1, k, n)) +
87 (Stencil::c22 * phi(i, j + 1, k - 1, n))) *
88 idx[2];
89 const amrex::Real phiz_jm1 =
90 ((Stencil::c20 * phi(i, j - 1, k + 1, n)) +
91 (Stencil::c21 * phi(i, j - 1, k, n)) +
92 (Stencil::c22 * phi(i, j - 1, k - 1, n))) *
93 idx[2];
94 const amrex::Real phiy =
95 ((Stencil::c10 * phi(i, j + 1, k, n)) +
96 (Stencil::c11 * phi(i, j, k, n)) +
97 (Stencil::c12 * phi(i, j - 1, k, n))) *
98 idx[1];
99 const amrex::Real phiy_ip1 =
100 ((Stencil::c10 * phi(i + 1, j + 1, k, n)) +
101 (Stencil::c11 * phi(i + 1, j, k, n)) +
102 (Stencil::c12 * phi(i + 1, j - 1, k, n))) *
103 idx[1];
104 const amrex::Real phiy_im1 =
105 ((Stencil::c10 * phi(i - 1, j + 1, k, n)) +
106 (Stencil::c11 * phi(i - 1, j, k, n)) +
107 (Stencil::c12 * phi(i - 1, j - 1, k, n))) *
108 idx[1];
109 const amrex::Real phiyz =
110 ((Stencil::c10 * phiz_jp1) + (Stencil::c11 * phiz) +
111 (Stencil::c12 * phiz_jm1)) *
112 idx[1];
113 const amrex::Real phix =
114 ((Stencil::c00 * phi(i + 1, j, k, n)) +
115 (Stencil::c01 * phi(i, j, k, n)) +
116 (Stencil::c02 * phi(i - 1, j, k, n))) *
117 idx[0];
118 const amrex::Real phixy =
119 ((Stencil::c00 * phiy_ip1) + (Stencil::c01 * phiy) +
120 (Stencil::c02 * phiy_im1)) *
121 idx[0];
122 const amrex::Real phixz =
123 ((Stencil::c00 * phiz_ip1) + (Stencil::c01 * phiz) +
124 (Stencil::c02 * phiz_im1)) *
125 idx[0];
126 curphi(i, j, k, n) =
127 -((phix * phix * phiyy) -
128 (2.0_rt * phix * phiy * phixy) +
129 (phiy * phiy * phixx) + (phix * phix * phizz) -
130 (2.0_rt * phix * phiz * phixz) +
131 (phiz * phiz * phixx) + (phiy * phiy * phizz) -
132 (2.0_rt * phiy * phiz * phiyz) +
133 (phiz * phiz * phiyy)) /
134 std::pow(
135 (phix * phix) + (phiy * phiy) + (phiz * phiz),
136 1.5_rt);
137 }
138 });
139 }
140
141 FTypeOut& m_curphi;
142 const FTypeIn& m_phi;
143};
144
148template <typename FTypeIn, typename FTypeOut>
149void curvature(FTypeOut& curphi, const FTypeIn& phi)
150{
151 BL_PROFILE("kynema-sgf::fvm::curvature");
152 Curvature<FTypeIn, FTypeOut> cur(curphi, phi);
153 impl::apply(cur, phi);
154}
155
159template <typename FType>
160std::unique_ptr<ScratchField> curvature(const FType& phi)
161{
162 const std::string gname = phi.name() + "_curvature";
163 auto curphi = phi.repo().create_scratch_field(gname, phi.num_comp());
164 curvature(*curphi, phi);
165 return curphi;
166}
167
168} // namespace kynema_sgf::fvm
169
170#endif /* CURVATURE_H */
fvm_utils.H
kynema_sgf::fvm::curvature
void curvature(FTypeOut &curphi, const FTypeIn &phi)
Definition curvature.H:149
kynema_sgf::fvm::impl::apply
AMREX_INLINE void apply(const FvmOp &fvmop, const FType &fld)
Definition fvm_utils.H:20
kynema_sgf::fvm
Definition SchemeTraits.H:6
kynema_sgf::fvm::Curvature
Definition curvature.H:16
kynema_sgf::fvm::Curvature::apply
void apply(const int lev) const
Definition curvature.H:33
kynema_sgf::fvm::Curvature::Curvature
Curvature(FTypeOut &curphi, const FTypeIn &phi)
Definition curvature.H:17
kynema_sgf::fvm::Curvature::m_phi
const FTypeIn & m_phi
Definition curvature.H:142
kynema_sgf::fvm::Curvature::m_curphi
FTypeOut & m_curphi
Definition curvature.H:141
Generated by doxygen 1.14.0