Line data Source code
1 : // Initialize a perturbed interface using Fourier Modes
2 : ///
3 : // Notes:
4 : // 1. \todo Extend this class to allow for 3D perturbations, currently only 2D are allowed
5 : // 2. \todo Allow for cosine (or complex exponential) expansions rather than just sin.
6 : // 3. \note This is a **two grain only** initial condition.
7 : // 4. \note This replaces the depricated "perturbed_bar" initial condition from previous versions
8 : //
9 : // The interface is defined as the :math:`x=0` plane (2D), or the :math:`x=0,z=0` plane (3D).
10 : // The equation for the interface is given by
11 : // :math:`y(x,z) = \sum_{n\in \{n_1,\ldots,n_N\}} A_n \sin(n\pi x/L_x)`
12 : // where :math:`A_n` are the amplitudes (stored in #wave_amplitudes),
13 : // :math:`n_1,\ldots,n_N\subset\mathbb{Z}_+` are wave numbers (stored in #wave_numbers),
14 : // and :math:`L_x` is the length in the x direction (obtained using the #geom object).
15 : //
16 : // Grain 1 is defined as being above :math:`y(x,z)`, Grain 2 is defined as being below.
17 :
18 : #ifndef IC_PERTURBEDINTERFACE_H_
19 : #define IC_PERTURBEDINTERFACE_H_
20 :
21 : #include "IC/IC.H"
22 : #include "Util/Util.H"
23 : #include "IO/ParmParse.H"
24 :
25 : namespace IC
26 : {
27 : class PerturbedInterface : public IC<Set::Scalar>
28 : {
29 : public:
30 : static constexpr const char* name = "perturbedinterface";
31 :
32 : enum Mollifier {Dirac, Gaussian};
33 0 : PerturbedInterface (amrex::Vector<amrex::Geometry> &_geom) :
34 0 : IC(_geom)
35 0 : { }
36 2 : PerturbedInterface(amrex::Vector<amrex::Geometry> &_geom,IO::ParmParse &pp, std::string name) : IC(_geom)
37 8 : {pp_queryclass(name,*this);}
38 : PerturbedInterface(amrex::Vector<amrex::Geometry> &_geom,IO::ParmParse &pp) : IC(_geom)
39 : {pp_queryclass(*this);}
40 :
41 20 : void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
42 : {
43 20 : Set::Scalar AMREX_D_DECL(L1 = geom[lev].ProbHi()[0] - geom[lev].ProbLo()[0],
44 : L2 = geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1],
45 : L3 = geom[lev].ProbHi()[2] - geom[lev].ProbLo()[2]);
46 :
47 1204 : for (amrex::MFIter mfi(*a_field[lev],true); mfi.isValid(); ++mfi)
48 : {
49 1184 : amrex::Box bx = mfi.tilebox();
50 1184 : bx.grow(a_field[lev]->nGrow());
51 1184 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
52 1184 : amrex::IndexType type = a_field[lev]->ixType();
53 :
54 1184 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
55 170496 : amrex::IntVect m(AMREX_D_DECL(i,j,k));
56 :
57 170496 : Set::Vector x;
58 : // NODE
59 340992 : if (type == amrex::IndexType::TheNodeType())
60 : {
61 0 : AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * geom[lev].CellSize()[0];,
62 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * geom[lev].CellSize()[1];,
63 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * geom[lev].CellSize()[2];);
64 : }
65 340992 : else if (type == amrex::IndexType::TheCellType())
66 : {
67 170496 : AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
68 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
69 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
70 : }
71 :
72 170496 : Set::Scalar bdry = 0.;
73 : #if AMREX_SPACEDIM == 2
74 170496 : Set::Scalar s1=NAN, l1=NAN;
75 170496 : if (normal==Direction::X) {s1 = x(1); l1 = L2;}
76 170496 : if (normal==Direction::Y) {s1 = x(0); l1 = L1;}
77 340992 : for (int n = 0; n < wave_numbers.size(); n++)
78 170496 : bdry += wave_amplitudes[n]
79 170496 : * ( fabs(std::cos(phis[n]))*std::cos(wave_numbers[n].real()*Set::Constant::Pi*s1 / l1) +
80 170496 : fabs(std::sin(phis[n]))*std::sin(wave_numbers[n].imag()*Set::Constant::Pi*s1 / l1));
81 : #elif AMREX_SPACEDIM == 3
82 0 : Set::Scalar s1=NAN, s2=NAN, l1=NAN, l2=NAN;
83 0 : if (normal==Direction::X) {s1 = x(1); s2 = x(2); l1 = L2; l2 = L3;}
84 0 : if (normal==Direction::Y) {s1 = x(2); s2 = x(0); l1 = L3; l2 = L1;}
85 0 : if (normal==Direction::Z) {s1 = x(0); s2 = x(0); l1 = L1; l2 = L2;}
86 0 : for (int n = 0; n < wave_numbers.size(); n++)
87 0 : bdry += wave_amplitudes[n]
88 0 : * ( fabs(std::cos(phis[n]))*std::cos(wave_numbers[n].real()*Set::Constant::Pi*s1 / l1) +
89 0 : fabs(std::sin(phis[n]))*std::sin(wave_numbers[n].imag()*Set::Constant::Pi*s1 / l1))
90 0 : * ( fabs(std::cos(phis[n]))*std::cos(wave_numbers[n].real()*Set::Constant::Pi*s2 / l2) +
91 0 : fabs(std::sin(phis[n]))*std::sin(wave_numbers[n].imag()*Set::Constant::Pi*s2 / l2))
92 : ;
93 : #endif
94 170496 : if (mol == Mollifier::Dirac)
95 : {
96 : // Util::Message(INFO);
97 0 : if ((normal == Direction::X && x(0) < bdry + offset)||
98 255744 : (normal == Direction::Y && x(1) < bdry + offset)||
99 85248 : (normal == Direction::Z && x(2) < bdry + offset))
100 : {
101 85248 : field(i,j,k,reverse) = 1.;
102 170496 : field(i,j,k,1-reverse) = 0.;
103 : }
104 : else
105 : {
106 85248 : field(i,j,k,reverse) = 0.;
107 170496 : field(i,j,k,1-reverse) = 1.;
108 : }
109 : }
110 : else
111 : {
112 0 : Set::Scalar t = 0.0;
113 0 : if (normal == Direction::X) t = x(0) - bdry - offset;
114 0 : else if (normal == Direction::Y) t = x(1) - bdry - offset;
115 0 : else if (normal == Direction::Z) t = x(2) - bdry - offset;
116 :
117 0 : Set::Scalar value = 0.5 + 0.5*std::erf(t/eps);
118 : // Util::Message(INFO, "value = ", value);
119 0 : field(i,j,k,reverse) = value;
120 0 : field(i,j,k,1-reverse) = 1. - value;
121 :
122 0 : if (field(i,j,k,0) < 0.0) field(i,j,k,reverse) = 0.0;
123 0 : if (field(i,j,k,0) > 1.0) field(i,j,k,reverse) = 1.0;
124 0 : if (field(i,j,k,1) < 0.0) field(i,j,k,1-reverse) = 0.0;
125 0 : if (field(i,j,k,1) > 1.0) field(i,j,k,1-reverse) = 1.0;
126 : }
127 170496 : });
128 20 : }
129 :
130 20 : };
131 :
132 : private:
133 : enum Direction {X,Y,Z};
134 : Direction normal = Direction::Y;
135 : Set::Scalar offset = 0.0;
136 : amrex::Vector<std::complex<int> > wave_numbers; ///< Store mode amplitudes \f$A_n\f$
137 : amrex::Vector<Set::Scalar> wave_amplitudes;
138 : std::vector<Set::Scalar> phis;
139 : Mollifier mol = Mollifier::Gaussian;
140 : Set::Scalar eps;
141 : int reverse = 0;
142 :
143 : public:
144 2 : static void Parse(PerturbedInterface & value, IO::ParmParse & pp)
145 : {
146 2 : std::vector<std::string> wave_numbers_str;
147 10 : pp_queryarr("wave_numbers",wave_numbers_str); // Wave numbers
148 4 : for (unsigned int i = 0; i<wave_numbers_str.size(); ++i)
149 : {
150 2 : value.wave_numbers.push_back(Util::String::Parse<std::complex<int> >(wave_numbers_str[i]));
151 2 : value.phis.push_back(std::atan2(value.wave_numbers[i].imag(),value.wave_numbers[i].real()));
152 : }
153 10 : pp_queryarr("wave_amplitudes",value.wave_amplitudes); // Wave amplitudes
154 :
155 2 : std::string normal_str;
156 2 : pp_query("normal",normal_str); // Which axis is normal to the interface (x,y,z)
157 2 : if (normal_str == "x") value.normal = Direction::X;
158 2 : if (normal_str == "y") value.normal = Direction::Y;
159 2 : if (normal_str == "z") value.normal = Direction::Z;
160 2 : std::string offset_str;
161 2 : pp_query("offset",value.offset); // Interface offset from origin
162 :
163 2 : pp_query("reverse",value.reverse); // If true, flip the interface (default:false)
164 14 : Util::Assert(INFO,TEST(value.reverse==0 || value.reverse==1));
165 :
166 2 : if (value.wave_numbers.size() != value.wave_amplitudes.size())
167 0 : Util::Abort(INFO, "Number of wave numbers and amplitudes must match");
168 :
169 2 : std::string mollifier = "dirac";
170 2 : pp_query("mollifier", mollifier); // Mollifier (options: dirac, [gaussian])
171 :
172 2 : if(mollifier == "dirac")
173 : {
174 2 : value.mol = Mollifier::Dirac;
175 4 : if (pp.contains("eps")) Util::Warning(INFO,"eps defined but not needed for dirac mollifier");
176 : }
177 : else
178 : {
179 0 : value.mol = Mollifier::Gaussian;
180 0 : pp_query("eps",value.eps); // Magnitude of mollifier
181 0 : Util::Message(INFO, "eps = ", value.eps);
182 : }
183 2 : }
184 :
185 : };
186 : }
187 : #endif
|
