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