Line data Source code
1 : //
2 : // This is an old-fashioned IC.
3 : // It is kept because it is still used with the Trig regression test, but
4 : // it is recommended to use IC::Expression instead.
5 : //
6 :
7 : #ifndef IC_TRIG_H_
8 : #define IC_TRIG_H_
9 :
10 : #include <complex>
11 :
12 : #include "IO/ParmParse.H"
13 : #include "AMReX_Vector.H"
14 : #include "IC/IC.H"
15 : #include "Util/Util.H"
16 : #include "Set/Set.H"
17 :
18 : namespace IC
19 : {
20 : /// \brief Initialize using a trigonometric series
21 : class Trig : public IC
22 : {
23 : public:
24 : static constexpr const char* name = "trig";
25 : Trig (amrex::Vector<amrex::Geometry> &_geom,
26 : Set::Scalar _alpha = 1.0,
27 : AMREX_D_DECL(std::complex<int> _n1 = 0,
28 : std::complex<int> _n2 = 0,
29 : std::complex<int> _n3 = 0),
30 : int _dim = AMREX_SPACEDIM) :
31 : IC(_geom)
32 : {
33 : Define(_alpha,AMREX_D_DECL(_n1,_n2,_n3),_dim);
34 : }
35 4 : Trig(amrex::Vector<amrex::Geometry> &_geom,IO::ParmParse &pp, std::string name) : IC(_geom)
36 4 : {pp_queryclass(name,*this);}
37 :
38 :
39 : void Define(Set::Scalar _alpha = 1.0,
40 : AMREX_D_DECL(std::complex<int> _n1 = 0,
41 : std::complex<int> _n2 = 0,
42 : std::complex<int> _n3 = 0),
43 : int _dim = AMREX_SPACEDIM)
44 : {
45 : alpha = _alpha;
46 : AMREX_D_TERM(n1 = _n1;, n2 = _n2;, n3 = _n3;);
47 : dim = _dim;
48 : AMREX_D_DECL(phi1 = std::atan2(n1.imag(),n1.real()),
49 : phi2 = std::atan2(n2.imag(),n2.real()),
50 : phi3 = std::atan2(n3.imag(),n3.real()));
51 : }
52 :
53 :
54 0 : void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
55 : {
56 0 : bool cellcentered = (a_field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
57 :
58 0 : const Set::Scalar AMREX_D_DECL(L1 = geom[lev].ProbHi()[0] - geom[lev].ProbLo()[0],
59 : L2 = geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1],
60 : L3 = geom[lev].ProbHi()[2] - geom[lev].ProbLo()[2]);
61 :
62 :
63 0 : for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
64 : {
65 0 : amrex::Box bx;
66 0 : if (cellcentered) bx = mfi.growntilebox();
67 0 : else bx = mfi.grownnodaltilebox();
68 :
69 0 : amrex::Array4<Set::Scalar> const &field = a_field[lev]->array(mfi);
70 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
71 : {
72 : Set::Scalar AMREX_D_DECL(x1,x2,x3);
73 0 : if (cellcentered)
74 : {
75 0 : AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
76 : x2 = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
77 : x3 = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
78 : }
79 : else
80 : {
81 0 : AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + (amrex::Real)(i) * geom[lev].CellSize()[0];,
82 : x2 = geom[lev].ProbLo()[1] + (amrex::Real)(j) * geom[lev].CellSize()[1];,
83 : x3 = geom[lev].ProbLo()[2] + (amrex::Real)(k) * geom[lev].CellSize()[2];);
84 : }
85 :
86 0 : Set::Scalar trigfn = 1.0;
87 : #if AMREX_SPACEDIM > 0
88 0 : if (dim > 0)
89 0 : trigfn *= (fabs(std::cos(phi1))*std::cos(n1.real()*Set::Constant::Pi*x1 / L1) +
90 0 : fabs(std::sin(phi1))*std::sin(n1.imag()*Set::Constant::Pi*x1 / L1));
91 : #endif
92 : #if AMREX_SPACEDIM > 1
93 0 : if (dim > 1)
94 0 : trigfn *= (fabs(std::cos(phi2))*std::cos(n2.real()*Set::Constant::Pi*x2 / L2) +
95 0 : fabs(std::sin(phi2))*std::sin(n2.imag()*Set::Constant::Pi*x2 / L2));
96 : #endif
97 : #if AMREX_SPACEDIM > 2
98 0 : if (dim > 2)
99 0 : trigfn *= (fabs(std::cos(phi3))*std::cos(n3.real()*Set::Constant::Pi*x3 / L3) +
100 0 : fabs(std::sin(phi3))*std::sin(n3.imag()*Set::Constant::Pi*x3 / L3));
101 : #endif
102 :
103 0 : field(i,j,k,comp) += alpha * trigfn;
104 :
105 0 : });
106 : }
107 0 : }
108 14 : void Add(const int &lev, Set::Field<Set::Vector> &a_field, Set::Scalar)
109 : {
110 14 : bool cellcentered = (a_field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
111 :
112 14 : const Set::Scalar AMREX_D_DECL(L1 = geom[lev].ProbHi()[0] - geom[lev].ProbLo()[0],
113 : L2 = geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1],
114 : L3 = geom[lev].ProbHi()[2] - geom[lev].ProbLo()[2]);
115 :
116 :
117 28 : for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
118 : {
119 14 : amrex::Box bx;
120 14 : if (cellcentered) bx = mfi.growntilebox();
121 14 : else bx = mfi.grownnodaltilebox();
122 :
123 14 : amrex::Array4<Set::Vector> const &field = a_field[lev]->array(mfi);
124 14 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
125 : {
126 : Set::Scalar AMREX_D_DECL(x1,x2,x3);
127 19166 : if (cellcentered)
128 : {
129 0 : AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
130 : x2 = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
131 : x3 = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
132 : }
133 : else
134 : {
135 19166 : AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + (amrex::Real)(i) * geom[lev].CellSize()[0];,
136 : x2 = geom[lev].ProbLo()[1] + (amrex::Real)(j) * geom[lev].CellSize()[1];,
137 : x3 = geom[lev].ProbLo()[2] + (amrex::Real)(k) * geom[lev].CellSize()[2];);
138 : }
139 :
140 19166 : Set::Scalar trigfn = 1.0;
141 : #if AMREX_SPACEDIM > 0
142 19166 : if (dim > 0)
143 19166 : trigfn *= (fabs(std::cos(phi1))*std::cos(n1.real()*Set::Constant::Pi*x1 / L1) +
144 19166 : fabs(std::sin(phi1))*std::sin(n1.imag()*Set::Constant::Pi*x1 / L1));
145 : #endif
146 : #if AMREX_SPACEDIM > 1
147 19166 : if (dim > 1)
148 19166 : trigfn *= (fabs(std::cos(phi2))*std::cos(n2.real()*Set::Constant::Pi*x2 / L2) +
149 19166 : fabs(std::sin(phi2))*std::sin(n2.imag()*Set::Constant::Pi*x2 / L2));
150 : #endif
151 : #if AMREX_SPACEDIM > 2
152 0 : if (dim > 2)
153 0 : trigfn *= (fabs(std::cos(phi3))*std::cos(n3.real()*Set::Constant::Pi*x3 / L3) +
154 0 : fabs(std::sin(phi3))*std::sin(n3.imag()*Set::Constant::Pi*x3 / L3));
155 : #endif
156 :
157 38332 : field(i,j,k)(comp) += alpha * trigfn;
158 :
159 19166 : });
160 : }
161 14 : }
162 : using IC::Add;
163 :
164 : public:
165 4 : static void Parse(Trig & value, IO::ParmParse & pp)
166 : {
167 8 : std::vector<int> n_real(AMREX_SPACEDIM,0.0);
168 4 : pp_queryarr("nr",n_real); // Number of real (cosin) waves
169 8 : std::vector<int> n_imag(AMREX_SPACEDIM,0.0);
170 4 : pp_queryarr("ni",n_imag); // Number of imaginary (sin) waves
171 :
172 4 : AMREX_D_TERM(
173 : value.n1 = std::complex<int>(n_real[0],n_imag[0]);,
174 : value.n2 = std::complex<int>(n_real[1],n_imag[1]);,
175 : value.n3 = std::complex<int>(n_real[2],n_imag[2]););
176 :
177 4 : AMREX_D_DECL(
178 : value.phi1 = std::atan2(value.n1.imag(),value.n1.real()),
179 : value.phi2 = std::atan2(value.n2.imag(),value.n2.real()),
180 : value.phi3 = std::atan2(value.n3.imag(),value.n3.real()));
181 :
182 4 : pp_query("dim",value.dim); // Spatial dimension
183 4 : pp_query("alpha",value.alpha); // Multiplier
184 :
185 4 : }
186 :
187 : private:
188 : int dim = AMREX_SPACEDIM;
189 : Set::Scalar alpha = 1.0;
190 : std::complex<int> AMREX_D_DECL(n1, n2, n3);
191 : Set::Scalar AMREX_D_DECL(phi1=0.0,phi2=0.0,phi3=0.0);
192 : };
193 : }
194 : #endif
|
