Line data Source code
1 : // This IC initializes a single-component fab using the formula
2 : //
3 : // .. math::
4 : //
5 : // \phi = \begin{cases}
6 : // \alpha_{in} & (\mathbf{x}-\mathbf{x}_0)^T\mathbf{A}(\mathbf{x}-\mathbf{x}_0) \le r \\
7 : // \alpha_{out} & \text{else}
8 : // \end{cases}
9 : //
10 : // The boundary can be mollified using an error function with parameter :math:`\varepsilon`
11 : //
12 : // .. WARNING::
13 : //
14 : // This IC is redundant with :ref:`IC::Ellipse` and :ref:`IC::Expression` and will probably be consolidated
15 : //
16 :
17 : #ifndef IC_ELLIPSOID_H_
18 : #define IC_ELLIPSOID_H_
19 :
20 : #include "IC/IC.H"
21 : #include "Util/Util.H"
22 : #include "IO/ParmParse.H"
23 :
24 : /// \class Ellipsoid
25 : /// \brief Initialize an ellipsoidal inclusion
26 : namespace IC
27 : {
28 : class Ellipsoid : public IC
29 : {
30 : public:
31 : enum Mollifier {Dirac, Gaussian};
32 :
33 0 : Ellipsoid (amrex::Vector<amrex::Geometry> &_geom) : IC(_geom) {}
34 :
35 0 : void Add(const int &lev,Set::Field<Set::Scalar> &a_field, Set::Scalar)
36 : {
37 0 : for (amrex::MFIter mfi(*a_field[lev],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
38 : {
39 0 : amrex::Box bx = mfi.tilebox();
40 0 : bx.grow(a_field[lev]->nGrow());
41 0 : amrex::IndexType type = a_field[lev]->ixType();
42 :
43 0 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
44 0 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
45 :
46 0 : Set::Vector x = Set::Position(i,j,k,geom[lev],type);
47 :
48 0 : Set::Scalar min_value = field(i,j,k);
49 0 : for (int i=0 ; i < center.size(); i++)
50 : {
51 0 : Set::Scalar value = 0.;
52 0 : Set::Scalar norm = (A[i]*(x-center[i])).lpNorm<2>();
53 0 : value = 0.5 + 0.5*std::erf(((x-center[i]).transpose() * A[i] * (x-center[i]) - 1.0) / eps[i] / norm);
54 0 : value = in_value + value*(out_value - in_value);
55 0 : min_value = value < min_value ? value : min_value;
56 : }
57 0 : field(i,j,k) = min_value;
58 0 : if (field(i,j,k) < std::min(in_value,out_value)) field(i,j,k) = std::min(in_value,out_value);
59 0 : if (field(i,j,k) > std::max(in_value,out_value)) field(i,j,k) = std::max(in_value,out_value);
60 0 : });
61 : }
62 0 : a_field[lev]->FillBoundary();
63 0 : };
64 :
65 : private:
66 : Set::Scalar in_value = 0.0, out_value = 1.0;
67 : amrex::Vector<Set::Vector> center;
68 : amrex::Vector<Set::Scalar> eps;
69 : amrex::Vector<Set::Matrix> A;
70 : Mollifier moll = Mollifier::Dirac;
71 :
72 : public:
73 0 : static void Parse(Ellipsoid & value, IO::ParmParse & pp)
74 : {
75 0 : amrex::Vector<Set::Scalar> center;
76 0 : if(pp.contains("center")) pp_queryarr("center",center); // Center of the ellipse :math:`\mathbf{x}_0`
77 :
78 0 : if(center.size() < AMREX_SPACEDIM) value.center.push_back(Set::Vector::Zero());
79 : else
80 : {
81 0 : for (int i = 0; i<center.size(); i+=AMREX_SPACEDIM)
82 0 : value.center.push_back(Set::Vector(AMREX_D_DECL(center[i],center[i+1],center[i+2])));
83 : }
84 :
85 0 : amrex::Vector<Set::Scalar> A;
86 0 : if (pp.contains("A"))
87 : {
88 0 : pp_queryarr("A", A); // Matrix defining elipse radii and orientation
89 0 : if(A.size() < AMREX_SPACEDIM*AMREX_SPACEDIM) value.A.push_back(Set::Matrix::Identity());
90 : else
91 : {
92 0 : int i=0, j=0;
93 0 : Set::Matrix temp = Set::Matrix::Zero();
94 0 : for (i=0; i<A.size(); i+=AMREX_SPACEDIM)
95 : {
96 0 : for (j=0; j<AMREX_SPACEDIM; j++)
97 0 : temp(i,j) = A[i+j];
98 :
99 0 : if((i+j) % (AMREX_SPACEDIM*AMREX_SPACEDIM) == 0)
100 : {
101 0 : value.A.push_back(temp);
102 0 : temp = Set::Matrix::Zero();
103 : }
104 : }
105 : }
106 : }
107 0 : if (pp.contains("radius"))
108 : {
109 0 : amrex::Vector<Set::Scalar> a_radius;
110 0 : pp_queryarr("radius",a_radius); // "Vector of radii (use instead of A)"
111 :
112 0 : if (a_radius.size() < AMREX_SPACEDIM) value.A.push_back(Set::Matrix::Identity());
113 : else
114 : {
115 0 : for (int i = 0; i<a_radius.size(); i+=AMREX_SPACEDIM)
116 : {
117 0 : Set::Matrix temp = Set::Matrix::Zero();
118 0 : AMREX_D_TERM( temp(0,0) = 1./(a_radius[i]*a_radius[i]) ;, temp(1,1) = 1./(a_radius[i+1]*a_radius[i+1]) ;, temp(2,2) = 1./(a_radius[i+2]*a_radius[i+2]) ; );
119 0 : value.A.push_back(temp);
120 : }
121 : }
122 : }
123 :
124 0 : amrex::Vector<Set::Scalar> eps;
125 0 : pp_queryarr("eps", eps); // Mollifying value for erf
126 :
127 0 : if(eps.size() < 1) value.eps.push_back(1.e-5);
128 0 : for (int i =0; i < value.A.size(); i++)
129 : {
130 0 : if (eps[i] <= 0.0) eps[i] = 1.e-5;
131 0 : value.eps.push_back(eps[i]);
132 : }
133 : //if(value.eps <= 0.) value.eps = 1.e-5;
134 :
135 0 : pp_query_default("in_value", value.in_value, 0.0); // Value of field inside ellipse
136 0 : pp_query_default("out_value", value.out_value, 1.0); // Value of field outside ellipse
137 :
138 0 : std::string mollifier;
139 0 : pp_query("mollifier",mollifier); // Type of mollifier to use (options: dirac, [gaussian])
140 0 : if(mollifier == "Dirac" || mollifier == "dirac")
141 0 : value.moll = Mollifier::Dirac;
142 : else
143 0 : value.moll = Mollifier::Gaussian;
144 0 : }
145 : };
146 : }
147 : #endif
|
