LCOV - code coverage report
Current view: top level - src/IC - Ellipse.H (source / functions) Hit Total Coverage
Test: coverage_merged.info Lines: 47 108 43.5 %
Date: 2025-02-18 04:54:05 Functions: 5 5 100.0 %

          Line data    Source code
       1             : // If :code:`number_of_inclusions` is specified, then multiple ellipses are specified.
       2             : // In this case, each parameter must have number_of_inclusion*M values, where M is the
       3             : // number of values specified for the single ellipse case.
       4             : //
       5             : 
       6             : #ifndef IC_ELLIPSE_H_
       7             : #define IC_ELLIPSE_H_
       8             : 
       9             : #include "Set/Set.H"
      10             : #include "IC/IC.H"
      11             : #include "IO/ParmParse.H"
      12             : 
      13             : namespace IC
      14             : {
      15             : class Ellipse : public IC
      16             : {
      17             : public:
      18             :     static constexpr const char* name = "ellipse";
      19             : 
      20             :     enum Mollifier {Dirac, Gaussian};
      21             : 
      22           5 :     Ellipse (amrex::Vector<amrex::Geometry> &_geom) : IC(_geom) {}
      23           5 :     Ellipse (amrex::Vector<amrex::Geometry> &_geom, IO::ParmParse &pp, std::string name) : Ellipse(_geom)
      24           5 :     {pp_queryclass(name,*this);}
      25             :     
      26          52 :     void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
      27             :     {
      28          52 :         Set::Vector DX(geom[lev].CellSize());
      29          52 :         amrex::IndexType type = a_field[lev]->ixType();
      30          52 :         int ncomp = a_field[lev]->nComp();
      31             :             
      32         313 :         for (amrex::MFIter mfi(*a_field[lev],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
      33             :         {
      34         261 :             amrex::Box bx;
      35         522 :             if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
      36         522 :             if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
      37             :             
      38         261 :             amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
      39         261 :             amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
      40             : 
      41      128591 :                 Set::Vector x;
      42             :                 // NODE
      43      257182 :                 if (type == amrex::IndexType::TheNodeType())
      44             :                 {
      45       77231 :                     AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * geom[lev].CellSize()[0];,
      46             :                                 x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * geom[lev].CellSize()[1];,
      47             :                                 x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * geom[lev].CellSize()[2];);
      48             :                 }
      49      102720 :                 else if (type == amrex::IndexType::TheCellType())
      50             :                 {
      51       51360 :                     AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
      52             :                                 x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
      53             :                                 x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
      54             :                 }
      55             : 
      56      128591 :                 if(number_of_inclusions == 0)
      57             :                 {   
      58      128591 :                     Set::Scalar norm = (A[0]*(x-center[0])).lpNorm<2>();
      59      128591 :                     field(i,j,k,0) = 0.5 - 0.5*std::erf(((x-center[0]).transpose() * A[0] * (x-center[0]) - 1.0) / eps[0] / norm);
      60      231311 :                     if (invert) field(i,j,k,0) = 1.0 - field(i,j,k,0);
      61      283053 :                     if (ncomp > 1) field(i,j,k,1) = 1.0 - field(i,j,k,0);
      62             :                 }
      63             :                 else
      64             :                 {
      65             :                     // Base matrix is 0, rest of the inclusions are numbered from 1.
      66           0 :                     Set::Scalar value = 0.0;
      67           0 :                     for (int m = 0; m < number_of_inclusions; m++)
      68             :                     {
      69           0 :                         Set::Scalar norm = (A[m]*(x-center[m])).lpNorm<2>();
      70           0 :                         value += 0.5 - 0.5*std::erf(((x-center[m]).transpose() * A[m] * (x-center[m]) - 1.0) / eps[m] / norm);
      71             :                         // if (field(i,j,k,m+1) < 0.) field(i,j,k,m+1) = 0.;
      72             :                         // if (field(i,j,k,m+1) > 1.) field(i,j,k,m+1) = 1.;
      73             :                         // value += field(i,j,k,m+1);
      74             :                     }
      75           0 :                     field (i,j,k,invert) = 1.0 - value;
      76           0 :                     field (i,j,k,1-invert) = value;
      77           0 :                     if (field(i,j,k,invert) < 0.) field(i,j,k,invert) = 0.;
      78           0 :                     if (field(i,j,k,invert) > 1.) field(i,j,k,invert) = 1.;
      79           0 :                     if (field(i,j,k,1-invert) < 0.) field(i,j,k,1-invert) = 0.;
      80           0 :                     if (field(i,j,k,1-invert) > 1.) field(i,j,k,1-invert) = 1.;
      81             :                 }
      82             :                 
      83      128591 :             });
      84             :         }
      85          52 :         a_field[lev]->FillBoundary();
      86          52 :     }
      87             :     
      88             : private:
      89             :     int number_of_inclusions = 1;
      90             :     amrex::Vector<Set::Vector> center;
      91             :     amrex::Vector<Set::Matrix> A;
      92             :     amrex::Vector<Set::Scalar> eps;
      93             :     int invert = 0;
      94             : 
      95             : public:
      96           5 :     static void Parse(Ellipse & value, IO::ParmParse & pp)
      97             :     {
      98          10 :         amrex::Vector<Set::Scalar> x0;
      99           5 :         if(!(pp.contains("number_of_inclusions")))
     100             :         {
     101           5 :             value.center.resize(0);
     102           5 :             value.A.resize(0);
     103           5 :             value.eps.resize(0);
     104             :             
     105           5 :             value.number_of_inclusions = 0;
     106           5 :             pp_queryarr("x0",x0); // Coorinates of ellipse center
     107           5 :             value.center.push_back(Set::Vector(AMREX_D_DECL(x0[0],x0[1],x0[2])));
     108             : 
     109             :             Set::Scalar _eps;
     110           5 :             pp_query_default("eps",_eps,0.0); // Diffuse boundary thickness
     111           5 :             value.eps.push_back(_eps);
     112             : 
     113           5 :             Set::Matrix _A = Set::Matrix::Zero();
     114           5 :             if (pp.contains("A"))
     115             :             {
     116           0 :                 pp_queryarr("A",_A); // DxD square matrix defining an ellipse. 
     117           0 :                 value.A.push_back(_A);
     118             :             }
     119           5 :             else if (pp.contains("a"))
     120             :             {
     121           5 :                 Set::Matrix _A = Set::Matrix::Zero();
     122           5 :                 Set::Vector a = Set::Vector::Ones();
     123           5 :                 pp_queryarr("a",a); // If :code:`A` is not defined, then assume a sphere with radius :code:`a`
     124          15 :                 for (int d = 0; d < AMREX_SPACEDIM; d++) _A(d,d) = 1./a(d)/a(d);
     125           5 :                 value.A.push_back(_A);
     126             :             }
     127             :         }
     128             :         else
     129             :         {
     130           0 :             pp_query("number_of_inclusions",value.number_of_inclusions); // Number of ellipses
     131           0 :             if(value.number_of_inclusions < 1) Util::Abort(INFO, "number of inclusions have to be at least 1");
     132             :             
     133           0 :             value.center.resize(0);
     134           0 :             value.A.resize(0);
     135           0 :             value.eps.resize(0);
     136             : 
     137           0 :             if (pp.contains("center") && pp.contains("x0")) Util::Abort(INFO,"Cannot specify both center (depricated) and x0");
     138           0 :             pp_queryarr("center", x0); // center of the ellipse
     139           0 :             pp_queryarr("x0", x0); // center of the ellipse
     140             :             
     141           0 :             if(x0.size() != value.number_of_inclusions*AMREX_SPACEDIM){
     142           0 :                 Util::Message(INFO, value.number_of_inclusions*AMREX_SPACEDIM);
     143           0 :                 Util::Message(INFO, x0.size());
     144           0 :                 Util::Abort(INFO, "Need centers for all the inclusions");}
     145           0 :             for (int i = 0; i < x0.size(); i+= AMREX_SPACEDIM)
     146           0 :                 value.center.push_back(Set::Vector(AMREX_D_DECL(x0[i],x0[i+1],x0[i+2])));
     147             : 
     148           0 :             if (pp.contains("A"))
     149             :             {
     150           0 :                 amrex::Vector<Set::Scalar> _A;
     151           0 :                 pp_queryarr("A", _A); // either a vector containing ellipse radii, or a matrix defining the ellipse
     152           0 :                 if(_A.size() != value.number_of_inclusions*AMREX_SPACEDIM*AMREX_SPACEDIM && _A.size() != AMREX_SPACEDIM*AMREX_SPACEDIM)
     153           0 :                     Util::Abort(INFO, "Invalid value of A for ellipse initialization");
     154           0 :                 if(_A.size() ==  AMREX_SPACEDIM*AMREX_SPACEDIM)
     155             :                 {
     156           0 :                     Set::Matrix _A1 = Set::Matrix::Zero();
     157           0 :                     pp_queryarr("A",_A1); // Same
     158           0 :                     for (int i = 0; i < value.number_of_inclusions; i++)
     159           0 :                         value.A.push_back(_A1);
     160             :                 } 
     161             :                 else
     162             :                 {
     163           0 :                     Set::Matrix _A1 = Set::Matrix::Zero();
     164           0 :                     for (int i = 0; i < value.number_of_inclusions; i+= AMREX_SPACEDIM*AMREX_SPACEDIM)
     165             :                     {
     166           0 :                         AMREX_D_PICK(   _A1(0,0) = _A[i];
     167             :                                         ,
     168             :                                         _A1(0,0) = _A[i]; _A1(0,1) = _A[i+1];
     169             :                                         _A1(1,0) = _A[i+2]; _A1(1,1) = _A[i+3];
     170             :                                         ,
     171             :                                         _A1(0,0) = _A[i+0]; _A1(0,1) = _A[i+1]; _A1(0,2) = _A[i+2];
     172             :                                         _A1(1,0) = _A[i+3]; _A1(1,1) = _A[i+4]; _A1(1,2) = _A[i+5];
     173             :                                         _A1(2,0) = _A[i+6]; _A1(2,1) = _A[i+7]; _A1(2,2) = _A[i+8];
     174             :                                     );
     175           0 :                         value.A.push_back(_A1);
     176             :                     }
     177             :                 }
     178             :             }
     179           0 :             else if (pp.contains("radius"))
     180             :             {
     181           0 :                 amrex::Vector<Set::Scalar> _r;
     182           0 :                 pp_queryarr("radius", _r); // Array of radii [depricated]
     183             : 
     184           0 :                 if(_r.size() != value.number_of_inclusions*AMREX_SPACEDIM && _r.size() != AMREX_SPACEDIM)
     185           0 :                     Util::Abort(INFO, "Invalid value of radius for ellipse initialization");
     186           0 :                 if(_r.size() ==  AMREX_SPACEDIM)
     187             :                 {
     188           0 :                     Set::Matrix _A1 = Set::Matrix::Zero();
     189           0 :                     for (int i = 0; i< AMREX_SPACEDIM; i++) _A1(i,i) = 1.0/(_r[i]*_r[i]);
     190           0 :                     for (int i = 0; i < value.number_of_inclusions; i++)
     191           0 :                         value.A.push_back(_A1);
     192             :                 } 
     193             :                 else
     194             :                 {
     195           0 :                     for (int i = 0; i < value.number_of_inclusions; i++)
     196             :                     {
     197           0 :                         Set::Matrix _A1 = Set::Matrix::Zero();
     198           0 :                         for (int d = 0; d <  AMREX_SPACEDIM; d++) _A1(d,d) = 1.0/_r[i+d]/_r[i+d];
     199           0 :                         value.A.push_back(_A1);
     200             :                     }
     201             :                 }
     202             :             }
     203             : 
     204           0 :             amrex::Vector<Set::Scalar> _eps;
     205           0 :             pp_queryarr("eps",_eps); // Regularization for smooth boundary
     206           0 :             if(_eps.size() != 1 && _eps.size() != value.number_of_inclusions)
     207           0 :                 Util::Abort(INFO, "Incorrect eps. Check the number of values specified");
     208             :             
     209           0 :             if (_eps.size() == 1)
     210           0 :                 for (int i = 0; i < value.number_of_inclusions; i++) value.eps.push_back(_eps[0]);
     211             :             
     212             :             else
     213           0 :                 for (int i = 0; i < value.number_of_inclusions; i++) value.eps.push_back(_eps[i]);
     214             :         }
     215             :         // Flip the inside and the outside
     216           5 :         pp_query("invert",value.invert);
     217           5 :     }
     218             : };
     219             : }
     220             : #endif

Generated by: LCOV version 1.14