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