LCOV - code coverage report
Current view: top level - src/IC - Laminate.H (source / functions) Hit Total Coverage
Test: coverage_merged.info Lines: 62 93 66.7 %
Date: 2024-11-18 05:28:54 Functions: 4 5 80.0 %

          Line data    Source code
       1             : //
       2             : // Create a single laminate with specified orientation, thickness, and offset.
       3             : //
       4             : 
       5             : #ifndef IC_LAMINATE_H_
       6             : #define IC_LAMINATE_H_
       7             : 
       8             : #include "IC/IC.H"
       9             : #include "Util/Util.H"
      10             : #include "IO/ParmParse.H"
      11             : 
      12             : /// \class Laminate
      13             : /// \brief Initialize Laminates in a matrix
      14             : namespace IC
      15             : {
      16             : class Laminate: public IC
      17             : {
      18             : public:
      19             :     static constexpr const char* name = "laminate";
      20             :     enum Mollifier { Dirac, Gaussian };
      21             : 
      22           0 :     Laminate(amrex::Vector<amrex::Geometry>& _geom): IC(_geom) {}
      23           3 :     Laminate(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name): IC(_geom)
      24             :     {        
      25           3 :         pp_queryclass(name, *this);
      26           3 :     }
      27             : 
      28          39 :     void Add(const int& lev, Set::Field<Set::Scalar>& a_field, Set::Scalar)
      29             :     {
      30         105 :         for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
      31             :         {
      32          66 :             amrex::Box bx = mfi.tilebox();
      33          66 :             bx.grow(a_field[lev]->nGrow());
      34          66 :             amrex::IndexType type = a_field[lev]->ixType();
      35             : 
      36          66 :             amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
      37          66 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
      38             : 
      39       73042 :                 Set::Vector x = Set::Position(i, j, k, geom[lev], type);
      40             : 
      41             :                 // The base matrix is always field(i,j,k,0)
      42             :                 // Inclusions start with field(i,j,k,1)
      43             : 
      44       73042 :                 if (!singlefab)
      45             :                 {
      46           0 :                     Set::Scalar value = 0.0;
      47           0 :                     for (int m = 0; m < number_of_inclusions; m++)
      48             :                     {
      49           0 :                         Set::Scalar t = std::abs((x - center[m]).transpose() * orientation[m]) - (thickness[m] / 2.0);
      50           0 :                         field(i, j, k, m + 1) = 0.5 - 0.5 * std::erf(t / eps[m]);
      51           0 :                         if (field(i, j, k, m + 1) < 0.) field(i, j, k, m + 1) = 0.;
      52           0 :                         if (field(i, j, k, m + 1) > 1.) field(i, j, k, m + 1) = 1.;
      53           0 :                         value += field(i, j, k, m + 1);
      54             :                     }
      55           0 :                     field(i, j, k, 0) = 1.0 - value;
      56           0 :                     if (field(i, j, k, 0) < 0.) field(i, j, k, 0) = 0.;
      57           0 :                     if (field(i, j, k, 0) > 1.) field(i, j, k, 0) = 1.;
      58             :                 }
      59             :                 else
      60             :                 {
      61       73042 :                     Set::Scalar t = std::abs((x - center[0]).transpose() * orientation[0]) - (thickness[0] / 2.0);
      62       73042 :                     field(i, j, k, 0) = 0.5 - 0.5 * std::erf(t / eps[0]);
      63      219126 :                     if (invert) field(i, j, k, 0) = 1.0 - field(i, j, k, 0);
      64             :                 }
      65       73042 :             });
      66             :         }
      67          39 :         a_field[lev]->FillBoundary();
      68          39 :     };
      69             : 
      70             : private:
      71             :     int number_of_inclusions = -1;
      72             :     amrex::Vector<Set::Vector> center;
      73             :     amrex::Vector<Set::Vector> orientation;
      74             :     amrex::Vector<Set::Vector> normal;
      75             :     amrex::Vector<Set::Scalar> eps;
      76             :     amrex::Vector<Set::Scalar> thickness;
      77             :     Mollifier moll = Mollifier::Gaussian;
      78             :     bool singlefab = false;
      79             :     bool invert = false;
      80             : 
      81             : public:
      82           3 :     static void Parse(Laminate& value, IO::ParmParse& pp)
      83             :     {
      84             :         // How many laminates (MUST be greater than or equal to 1).
      85           3 :         pp_query_default("number_of_inclusions", value.number_of_inclusions, 1);
      86             : 
      87           3 :         if (value.number_of_inclusions < 1)
      88           0 :             Util::Abort(INFO, "Number of inclusions must be at least 1. Aborting.");
      89             : 
      90           6 :         amrex::Vector<Set::Scalar> a_center;
      91             :         // (x,y,[z]) values for the center point of the laminate
      92           3 :         if (pp.contains("center")) pp_queryarr("center", a_center);
      93             : 
      94           3 :         if (a_center.size() != value.number_of_inclusions * AMREX_SPACEDIM) value.center.push_back(Set::Vector::Zero());
      95             :         else
      96             :         {
      97           0 :             for (int i = 0; i < a_center.size(); i += AMREX_SPACEDIM)
      98           0 :                 value.center.push_back(Set::Vector(AMREX_D_DECL(a_center[i], a_center[i + 1], a_center[i + 2])));
      99             :         }
     100             : 
     101           6 :         amrex::Vector<Set::Scalar> a_thickness;
     102             :         // thickness of the laminate
     103           3 :         if (pp.contains("thickness")) pp_queryarr("thickness", a_thickness);
     104             : 
     105           3 :         if (a_thickness.size() != value.number_of_inclusions && a_thickness.size() != 1)
     106           0 :             Util::Abort(INFO, "Thickness of each inclusion must be specified");
     107             : 
     108           3 :         if (a_thickness.size() == 1)
     109             :         {
     110           3 :             if (a_thickness[0] <= 0.0) Util::Abort(INFO, "Invalid value of inclusion thickness");
     111           6 :             for (int i = 0; i < value.number_of_inclusions; i++) value.thickness.push_back(a_thickness[0]);
     112             :         }
     113             : 
     114             :         else
     115             :         {
     116           0 :             for (int i = 0; i < value.number_of_inclusions; i++)
     117             :             {
     118           0 :                 if (a_thickness[i] <= 0.0) Util::Abort(INFO, "Invalid value of inclusion ", i + 1, " thickness");
     119           0 :                 value.thickness.push_back(a_thickness[i]);
     120             :             }
     121             :         }
     122             : 
     123           6 :         amrex::Vector<Set::Scalar> a_orientation;
     124             :         // Vector normal to the interface of the laminate
     125           3 :         pp_queryarr("orientation", a_orientation);
     126             : 
     127           3 :         if (a_orientation.size() != value.number_of_inclusions * AMREX_SPACEDIM && a_orientation.size() != AMREX_SPACEDIM)
     128           0 :             Util::Abort(INFO, "Orientation of each inclusion must be specified");
     129             : 
     130           3 :         if (a_orientation.size() == AMREX_SPACEDIM)
     131           6 :             for (int i = 0; i < value.number_of_inclusions; i++)
     132           3 :                 value.orientation.push_back(Set::Vector(AMREX_D_DECL(a_orientation[0], a_orientation[1], a_orientation[2])));
     133             : 
     134             :         else
     135           0 :             for (int i = 0; i < a_orientation.size(); i += AMREX_SPACEDIM)
     136           0 :                 value.orientation.push_back(Set::Vector(AMREX_D_DECL(a_orientation[i], a_orientation[i + 1], a_orientation[i + 2])));
     137             : 
     138           6 :         for (int i = 0; i < value.orientation.size(); i++)
     139           3 :             if (value.orientation[i].lpNorm<2>() <= 0.) value.orientation[i] = Set::Vector::Random();
     140             : 
     141           6 :         for (int i = 0; i < value.orientation.size(); i++)
     142           3 :             value.orientation[i] = value.orientation[i] / value.orientation[i].lpNorm<2>();
     143             : 
     144           3 :         value.normal.resize(value.number_of_inclusions);
     145           6 :         for (int i = 0; i < value.orientation.size(); i++)
     146             :         {
     147           3 :             value.normal[i] = Set::Vector::Zero();
     148           3 :             if (value.orientation[i](0) != 0.)
     149             :             {
     150           0 :                 AMREX_D_TERM(value.normal[i](0) = 1.;, value.normal[i](1) = 1.;, value.normal[i](2) = 1.;);
     151           0 :                 value.normal[i](0) = -(AMREX_D_TERM(0., +value.orientation[i](1), +value.orientation[i](2))) / value.orientation[i](0);
     152           0 :                 value.normal[i] = value.normal[i] / value.normal[i].lpNorm<2>();
     153             :             }
     154           3 :             else if (value.orientation[i](1) != 0.)
     155             :             {
     156           3 :                 AMREX_D_TERM(value.normal[i](0) = 1.;, value.normal[i](1) = 1.;, value.normal[i](2) = 1.;);
     157           3 :                 value.normal[i](1) = -(AMREX_D_TERM(value.orientation[i](0), +0.0, +value.orientation[i](2))) / value.orientation[i](1);
     158           3 :                 value.normal[i] = value.normal[i] / value.normal[i].lpNorm<2>();
     159             :             }
     160             :         }
     161             : 
     162           6 :         amrex::Vector<Set::Scalar> a_eps;
     163             :         // Diffuse thickness
     164           3 :         pp_queryarr("eps", a_eps);
     165             : 
     166           3 :         if (a_eps.size() < 1)
     167           0 :             for (int i = 0; i < value.number_of_inclusions; i++)
     168           0 :                 value.eps.push_back(1.e-5);
     169           3 :         if (a_eps.size() == 1)
     170             :         {
     171           3 :             if (a_eps[0] < 0.0)
     172             :             {
     173           0 :                 Util::Warning(INFO, "Invalid value of laminate.eps. Resetting to 1e-5");
     174           0 :                 a_eps[0] = 1.e-5;
     175             :             }
     176           6 :             for (int i = 0; i < value.number_of_inclusions; i++)
     177           3 :                 value.eps.push_back(a_eps[0]);
     178             :         }
     179             :         else
     180             :         {
     181           0 :             for (int i = 0; i < value.number_of_inclusions; i++)
     182           0 :                 value.eps.push_back(a_eps[i]);
     183             :         }
     184             : 
     185           6 :         std::string mollifier;
     186           3 :         pp_query("mollifier", mollifier); // Type of mollifer to use (options: dirac, [gaussian])
     187           3 :         if (mollifier == "Dirac" || mollifier == "dirac")
     188           0 :             value.moll = Mollifier::Dirac;
     189             :         else
     190           3 :             value.moll = Mollifier::Gaussian;
     191             : 
     192             :         // Switch to mode where only one component is used.
     193           3 :         pp_query("singlefab", value.singlefab);
     194             : 
     195             :         // Take the complement of the laminate
     196           3 :         pp_query("invert", value.invert);
     197           3 :     }
     198             : };
     199             : }
     200             : #endif

Generated by: LCOV version 1.14