Alamo
Trig.H
Go to the documentation of this file.
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  Trig(amrex::Vector<amrex::Geometry> &_geom,IO::ParmParse &pp, std::string name) : IC(_geom)
35  {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  void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
54  {
55  bool cellcentered = (a_field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
56 
57  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  for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
63  {
64  amrex::Box bx;
65  if (cellcentered) bx = mfi.growntilebox();
66  else bx = mfi.grownnodaltilebox();
67 
68  amrex::Array4<Set::Scalar> const &field = a_field[lev]->array(mfi);
69  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
70  {
71  Set::Scalar AMREX_D_DECL(x1,x2,x3);
72  if (cellcentered)
73  {
74  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  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  Set::Scalar trigfn = 1.0;
86  #if AMREX_SPACEDIM > 0
87  if (dim > 0)
88  trigfn *= (fabs(std::cos(phi1))*std::cos(n1.real()*Set::Constant::Pi*x1 / L1) +
89  fabs(std::sin(phi1))*std::sin(n1.imag()*Set::Constant::Pi*x1 / L1));
90  #endif
91  #if AMREX_SPACEDIM > 1
92  if (dim > 1)
93  trigfn *= (fabs(std::cos(phi2))*std::cos(n2.real()*Set::Constant::Pi*x2 / L2) +
94  fabs(std::sin(phi2))*std::sin(n2.imag()*Set::Constant::Pi*x2 / L2));
95  #endif
96  #if AMREX_SPACEDIM > 2
97  if (dim > 2)
98  trigfn *= (fabs(std::cos(phi3))*std::cos(n3.real()*Set::Constant::Pi*x3 / L3) +
99  fabs(std::sin(phi3))*std::sin(n3.imag()*Set::Constant::Pi*x3 / L3));
100  #endif
101 
102  field(i,j,k,comp) += alpha * trigfn;
103 
104  });
105  }
106  }
107  void Add(const int &lev, Set::Field<Set::Vector> &a_field, Set::Scalar)
108  {
109  bool cellcentered = (a_field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
110 
111  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  for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
117  {
118  amrex::Box bx;
119  if (cellcentered) bx = mfi.growntilebox();
120  else bx = mfi.grownnodaltilebox();
121 
122  amrex::Array4<Set::Vector> const &field = a_field[lev]->array(mfi);
123  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
124  {
125  Set::Scalar AMREX_D_DECL(x1,x2,x3);
126  if (cellcentered)
127  {
128  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  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  Set::Scalar trigfn = 1.0;
140  #if AMREX_SPACEDIM > 0
141  if (dim > 0)
142  trigfn *= (fabs(std::cos(phi1))*std::cos(n1.real()*Set::Constant::Pi*x1 / L1) +
143  fabs(std::sin(phi1))*std::sin(n1.imag()*Set::Constant::Pi*x1 / L1));
144  #endif
145  #if AMREX_SPACEDIM > 1
146  if (dim > 1)
147  trigfn *= (fabs(std::cos(phi2))*std::cos(n2.real()*Set::Constant::Pi*x2 / L2) +
148  fabs(std::sin(phi2))*std::sin(n2.imag()*Set::Constant::Pi*x2 / L2));
149  #endif
150  #if AMREX_SPACEDIM > 2
151  if (dim > 2)
152  trigfn *= (fabs(std::cos(phi3))*std::cos(n3.real()*Set::Constant::Pi*x3 / L3) +
153  fabs(std::sin(phi3))*std::sin(n3.imag()*Set::Constant::Pi*x3 / L3));
154  #endif
155 
156  field(i,j,k)(comp) += alpha * trigfn;
157 
158  });
159  }
160  }
161  using IC::Add;
162 
163 public:
164  static void Parse(Trig & value, IO::ParmParse & pp)
165  {
166  std::vector<int> n_real(AMREX_SPACEDIM,0.0);
167  pp_queryarr("nr",n_real); // Number of real (cosin) waves
168  std::vector<int> n_imag(AMREX_SPACEDIM,0.0);
169  pp_queryarr("ni",n_imag); // Number of imaginary (sin) waves
170 
171  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  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  pp_query("dim",value.dim); // Spatial dimension
182  pp_query("alpha",value.alpha); // Multiplier
183 
184  }
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
IC::IC::geom
amrex::Vector< amrex::Geometry > & geom
Definition: IC.H:45
IC::Trig
Initialize using a trigonometric series.
Definition: Trig.H:21
IC::IC::Add
virtual void Add(const int &lev, Set::Field< Set::Scalar > &field, Set::Scalar time)=0
Util.H
Set::Field< Set::Scalar >
Definition: Set.H:236
IC::Trig::Define
void Define(Set::Scalar _alpha=1.0, AMREX_D_DECL(std::complex< int > _n1=0, std::complex< int > _n2=0, std::complex< int > _n3=0), int _dim=AMREX_SPACEDIM)
Definition: Trig.H:38
ParmParse.H
pp_query
#define pp_query(...)
Definition: ParmParse.H:104
IC::Trig::dim
int dim
Definition: Trig.H:187
Set::Scalar
amrex::Real Scalar
Definition: Base.H:19
pp_queryclass
#define pp_queryclass(...)
Definition: ParmParse.H:105
IC::Trig::Parse
static void Parse(Trig &value, IO::ParmParse &pp)
Definition: Trig.H:164
IC::Trig::AMREX_D_DECL
std::complex< int > AMREX_D_DECL(n1, n2, n3)
IC::Trig::Trig
Trig(amrex::Vector< amrex::Geometry > &_geom, IO::ParmParse &pp, std::string name)
Definition: Trig.H:34
IC::IC::comp
int comp
Definition: IC.H:46
Set.H
Set::Constant::Pi
static const Set::Scalar Pi
Definition: Set.H:286
IC
Definition: Affine.H:18
IO::ParmParse
Definition: ParmParse.H:110
IC::Trig::alpha
Set::Scalar alpha
Definition: Trig.H:188
IC::Trig::Add
void Add(const int &lev, Set::Field< Set::Vector > &a_field, Set::Scalar)
Definition: Trig.H:107
IC::Trig::Trig
Trig(amrex::Vector< amrex::Geometry > &_geom, Set::Scalar _alpha=1.0, AMREX_D_DECL(std::complex< int > _n1=0, std::complex< int > _n2=0, std::complex< int > _n3=0), int _dim=AMREX_SPACEDIM)
Definition: Trig.H:24
IC::Trig::Add
void Add(const int &lev, Set::Field< Set::Scalar > &a_field, Set::Scalar)
Definition: Trig.H:53
IC.H
Set::Field< Set::Vector >
pp_queryarr
#define pp_queryarr(...)
Definition: ParmParse.H:103
Generated on Fri Aug 16 2024 16:30:29 for Alamo by   doxygen 1.8.17