MacroscopicProperties.H

Go to the documentation of this file.

/* Copyright 2020 Revathi Jambunathan
 *
 * This file is part of WarpX.
 *
 * License: BSD-3-Clause-LBNL
 */
 
#ifndef WARPX_MACROSCOPICPROPERTIES_H_
#define WARPX_MACROSCOPICPROPERTIES_H_
 
#include "MacroscopicProperties_fwd.H"
 
#include "Utils/WarpXConst.H"
 
#include <AMReX_Array.H>
#include <AMReX_Extension.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_MultiFab.H>
#include <AMReX_Parser.H>
#include <AMReX_REAL.H>
 
#include <memory>
#include <string>
 
 
class
MacroscopicProperties
{
public:
     MacroscopicProperties (); // constructor
     void ReadParameters ();
     void InitData ();
 
     amrex::MultiFab& getsigma_mf  () {return (*m_sigma_mf);}
     amrex::MultiFab * get_pointer_sigma () {return m_sigma_mf.get();}
     amrex::MultiFab& getepsilon_mf  () {return (*m_eps_mf);}
     amrex::MultiFab * get_pointer_eps () {return m_eps_mf.get();}
     amrex::MultiFab& getmu_mf  () {return (*m_mu_mf);}
     amrex::MultiFab * get_pointer_mu () {return m_mu_mf.get();}
 
     amrex::GpuArray<int, 3> macro_cr_ratio;
     void InitializeMacroMultiFabUsingParser (amrex::MultiFab *macro_mf,
                                              amrex::ParserExecutor<3> const& macro_parser,
                                              const int lev);
 
    void InitializeMacroMultiFabFromNumpy (amrex::MultiFab* macro_mf,
                                           const std::string& npy_filename,
                                           const int lev,
                                           const int npy_k_index,
                                           const amrex::Real npy_value);
 
    void InitializePECFromSigma (amrex::MultiFab* sigma_mf,
                                 amrex::MultiFab* PECx,
                                 amrex::MultiFab* PECy,
                                 const int npy_k_index);
 
     amrex::GpuArray<int, 3> sigma_IndexType;
     amrex::GpuArray<int, 3> epsilon_IndexType;
     amrex::GpuArray<int, 3> mu_IndexType;
     amrex::GpuArray<int, 3> Ex_IndexType;
     amrex::GpuArray<int, 3> Ey_IndexType;
     amrex::GpuArray<int, 3> Ez_IndexType;
     amrex::GpuArray<int, 3> Bx_IndexType;
     amrex::GpuArray<int, 3> By_IndexType;
     amrex::GpuArray<int, 3> Bz_IndexType;
 
     std::string m_sigma_s = "constant";
     std::string m_epsilon_s = "constant";
     std::string m_mu_s = "constant";
     amrex::Real m_sigma = 0.0;
     amrex::Real m_epsilon = PhysConst::ep0;
     amrex::Real m_mu = PhysConst::mu0;
     int m_npy_k_index = 0;
     std::string m_sigma_npy_filename;
     std::string m_epsilon_npy_filename;
     std::string m_mu_npy_filename;
     amrex::Real m_sigma_npy_value = 0.;
     amrex::Real m_epsilon_npy_value = 0.;
     amrex::Real m_mu_npy_value = 0.;
     std::unique_ptr<amrex::Parser> m_sigma_parser;
     std::unique_ptr<amrex::Parser> m_epsilon_parser;
     std::unique_ptr<amrex::Parser> m_mu_parser;
 
#ifdef WARPX_MAG_LLG
     amrex::GpuArray<int, 3> Hx_IndexType;
     amrex::GpuArray<int, 3> Hy_IndexType;
     amrex::GpuArray<int, 3> Hz_IndexType;
     amrex::GpuArray<int, 3> Mx_IndexType;
     amrex::GpuArray<int, 3> My_IndexType;
     amrex::GpuArray<int, 3> Mz_IndexType;
     amrex::GpuArray<amrex::Real, 3> mag_LLG_anisotropy_axis;
 
     amrex::MultiFab& getmag_Ms_mf        (int dir) {return (*m_mag_Ms_mf[dir]);}
     amrex::MultiFab * getmag_pointer_Ms (int dir) {return m_mag_Ms_mf[dir].get();}
     amrex::MultiFab& getmag_alpha_mf     (int dir) {return (*m_mag_alpha_mf[dir]);}
     amrex::MultiFab * getmag_pointer_alpha (int dir) {return m_mag_alpha_mf[dir].get();}
     amrex::MultiFab& getmag_gamma_mf     (int dir) {return (*m_mag_gamma_mf[dir]);}
     amrex::MultiFab * getmag_pointer_gamma (int dir) {return m_mag_gamma_mf[dir].get();}
     amrex::MultiFab& getmag_exchange_mf  (int dir) {return (*m_mag_exchange_mf[dir]);}
     amrex::MultiFab * getmag_pointer_exchange (int dir) {return m_mag_exchange_mf[dir].get();}
     amrex::MultiFab& getmag_anisotropy_mf(int dir) {return (*m_mag_anisotropy_mf[dir]);}
     amrex::MultiFab * getmag_pointer_anisotropy (int dir) {return m_mag_anisotropy_mf[dir].get();}
 
     amrex::Real getmag_normalized_error () {return m_mag_normalized_error;}
     int getmag_max_iter () {return m_mag_max_iter;}
     amrex::Real getmag_tol () {return m_mag_tol;}
 
    // interpolate the magnetic properties to B locations
     // magnetic properties are cell nodal
     // B locations are face centered
     // iv is an IntVect with a 1 in the face direction of interest, and 0 in the others
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real macro_avg_to_face (int i, int j, int k, amrex::IntVect iv, amrex::Array4<amrex::Real> const& macro_mag_prop){
         using namespace amrex;
         return ( 0.125_rt * ( macro_mag_prop(i        ,j        ,k        )
                             + macro_mag_prop(i-iv[0]+1,j        ,k        )
                             + macro_mag_prop(i        ,j-iv[1]+1,k        )
                             + macro_mag_prop(i        ,j        ,k-iv[2]+1)
                             + macro_mag_prop(i        ,j-iv[1]+1,k-iv[2]+1)
                             + macro_mag_prop(i-iv[0]+1,j        ,k-iv[2]+1)
                             + macro_mag_prop(i-iv[0]+1,j-iv[1]+1,k        )
                             + macro_mag_prop(i-iv[0]+1,j-iv[1]+1,k-iv[2]+1)
         ));
     }
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real face_avg_to_face (int i, int j, int k, int n,
                                           amrex::IntVect iv_in, amrex::IntVect iv_out,
                                           amrex::Array4<amrex::Real> const& Fieldcomp) {
         using namespace amrex;
         return ( 0.125_rt * ( Fieldcomp(i                   , j                   , k                   , n)
                             + Fieldcomp(i+iv_in[0]-iv_out[0], j                   , k                   , n)
                             + Fieldcomp(i                   , j+iv_in[1]-iv_out[1], k                   , n)
                             + Fieldcomp(i                   , j                   , k+iv_in[2]-iv_out[2], n)
                             + Fieldcomp(i+iv_in[0]-iv_out[0], j+iv_in[1]-iv_out[1], k                   , n)
                             + Fieldcomp(i+iv_in[0]-iv_out[0], j                   , k+iv_in[2]-iv_out[2], n)
                             + Fieldcomp(i                   , j+iv_in[1]-iv_out[1], k+iv_in[2]-iv_out[2], n)
                             + Fieldcomp(i+iv_in[0]-iv_out[0], j+iv_in[1]-iv_out[1], k+iv_in[2]-iv_out[2], n)
                             ));
     }
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real getH_Maxwell (int i, int j, int k, int n,
                                   amrex::IntVect iv_in, amrex::IntVect iv_out,
                                   amrex::Array4<amrex::Real> const& Bcomp, amrex::Array4<amrex::Real> const& Mcomp) {
         using namespace amrex;
         amrex::Real H_Maxwell = face_avg_to_face(i, j, k, 0, iv_in, iv_out, Bcomp)/PhysConst::mu0
                                 - Mcomp(i, j, k, n); //magnetic constitutive relation
         return H_Maxwell;
     }
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real updateM_field (int i, int j, int k, int n,
                                   amrex::Array4<amrex::Real> const& a, amrex::Array4<amrex::Real> const& b) {
         using namespace amrex;
         amrex::Real a_square = pow(a(i, j, k, 0), 2.0) + pow(a(i, j, k, 1), 2.0) + pow(a(i, j, k, 2), 2.0);
         amrex::Real a_dot_b =  a(i, j, k, 0) * b(i, j, k, 0) +
                                a(i, j, k, 1) * b(i, j, k, 1) +
                                a(i, j, k, 2) * b(i, j, k, 2);
         // Initialize to 0.
         amrex::Real M_field = 0._rt;
 
         if(n==0){
             amrex::Real a_cross_b_x = a(i, j, k, 1) * b(i, j, k, 2) -
                                       a(i, j, k, 2) * b(i, j, k, 1);
             M_field = ( b(i, j, k, 0) + a_dot_b * a(i, j, k, 0) - a_cross_b_x ) / ( 1.0 + a_square);
         }
         else if(n==1){
             amrex::Real a_cross_b_y = a(i, j, k, 2) * b(i, j, k, 0) -
                                       a(i, j, k, 0) * b(i, j, k, 2);
             M_field = ( b(i, j, k, 1) + a_dot_b * a(i, j, k, 1) - a_cross_b_y ) / ( 1.0 + a_square);
         }
         else if(n==2){
             amrex::Real a_cross_b_z = a(i, j, k, 0) * b(i, j, k, 1) -
                                       a(i, j, k, 1) * b(i, j, k, 0);
             M_field = ( b(i, j, k, 2) + a_dot_b * a(i, j, k, 2) - a_cross_b_z ) / ( 1.0 + a_square);
         }
         else{
             amrex::Abort("Wrong component n of the M_field");
         }
         return M_field;
     }
#endif //closes ifdef MAG_LLG
 
#ifdef WARPX_MAG_LLG
     amrex::Real m_mag_Ms;
     amrex::Real m_mag_alpha;
     amrex::Real m_mag_gamma;
     amrex::Real m_mag_exchange;
     amrex::Real m_mag_anisotropy;
 
     // If the magnitude of the magnetization deviates by more than this amount relative
     // to the user-defined Ms, abort.  Default 0.1.
     amrex::Real m_mag_normalized_error;
 
     // maximum iteration count for the second-order time advancement scheme of M field, default 100
     int m_mag_max_iter;
 
     // the relative tolerance for the second-order time advancement scheme of M field, default 0.0001
     amrex::Real m_mag_tol;
 
     std::array<std::unique_ptr<amrex::MultiFab>, 3> m_mag_Ms_mf;
     std::array<std::unique_ptr<amrex::MultiFab>, 3> m_mag_alpha_mf;
     std::array<std::unique_ptr<amrex::MultiFab>, 3> m_mag_gamma_mf;
     std::array<std::unique_ptr<amrex::MultiFab>, 3> m_mag_exchange_mf;
     std::array<std::unique_ptr<amrex::MultiFab>, 3> m_mag_anisotropy_mf;
 
     // these store the type of initialization, e.g., "constant", "parse_X_function", etc.
     std::string m_mag_Ms_s;
     std::string m_mag_alpha_s;
     std::string m_mag_gamma_s;
     std::string m_mag_exchange_s;
     std::string m_mag_anisotropy_s;
 
     // these store the parsed expression for the material properties
     std::string m_str_mag_Ms_function;
     std::string m_str_mag_alpha_function;
     std::string m_str_mag_gamma_function;
     std::string m_str_mag_exchange_function;
     std::string m_str_mag_anisotropy_function;
 
     std::unique_ptr<amrex::Parser> m_mag_Ms_parser;
     std::unique_ptr<amrex::Parser> m_mag_alpha_parser;
     std::unique_ptr<amrex::Parser> m_mag_gamma_parser;
     std::unique_ptr<amrex::Parser> m_mag_exchange_parser;
     std::unique_ptr<amrex::Parser> m_mag_anisotropy_parser;
#endif
 
private:
 
     std::unique_ptr<amrex::MultiFab> m_sigma_mf;
     std::unique_ptr<amrex::MultiFab> m_eps_mf;
     std::unique_ptr<amrex::MultiFab> m_mu_mf;
 
 
     std::string m_str_sigma_function;
     std::string m_str_epsilon_function;
     std::string m_str_mu_function;
 
};
 
struct LaxWendroffAlgo {
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real alpha (amrex::Real const sigma,
                               amrex::Real const epsilon,
                               amrex::Real dt) {
         using namespace amrex;
         amrex::Real fac1 = 0.5_rt * sigma * dt / epsilon;
         amrex::Real alpha = (1._rt - fac1)/(1._rt + fac1);
         return alpha;
     }
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real beta (amrex::Real const sigma,
                              amrex::Real const epsilon,
                              amrex::Real dt) {
         using namespace amrex;
         amrex::Real fac1 = 0.5_rt * sigma * dt / epsilon;
         amrex::Real beta = dt / ( epsilon * (1._rt + fac1) );
         return beta;
     }
 
};
 
struct BackwardEulerAlgo {
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real alpha (amrex::Real const sigma,
                               amrex::Real const epsilon,
                               amrex::Real dt) {
         using namespace amrex;
         amrex::Real fac1 = sigma * dt / epsilon;
         amrex::Real alpha = (1._rt)/(1._rt + fac1);
         return alpha;
     }
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real beta (amrex::Real const sigma,
                              amrex::Real const epsilon,
                              amrex::Real dt) {
         using namespace amrex;
         amrex::Real fac1 = sigma * dt / epsilon;
         amrex::Real beta = dt / ( epsilon * (1._rt + fac1) );
         return beta;
     }
 
};
 
#endif // WARPX_MACROSCOPIC_PROPERTIES_H_