initialization and diagnostics for polarization

Source/Diagnostics/Diagnostics.cpp

@@ -173,6 +173,24 @@ Diagnostics::BaseReadParameters ()
 
 #endif
 
+#ifdef WARPX_FERROE
+    // polarization can be written to file only if WarpX::em_solver_medium == MediumForEM::Macroscopic
+    if (utils::algorithms::is_in(m_varnames, "px") || utils::algorithms::is_in(m_varnames, "py") || utils::algorithms::is_in(m_varnames, "pz"))
+    {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            WarpX::em_solver_medium == MediumForEM::Macroscopic,
+            "polarization in plotfiles only works with algo.em_solver_medium=macroscopic");
+    }
+
+    // ferroelectric can be written to file only if WarpX::em_solver_medium == MediumForEM::Macroscopic
+    if (utils::algorithms::is_in(m_varnames, "ferroelectric"))
+    {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            WarpX::em_solver_medium == MediumForEM::Macroscopic,
+            "ferroelectric in plotfiles only works with algo.em_solver_medium=macroscopic");
+    }
+#endif
+
     // If user requests to plot proc_number for a serial run,
     // delete proc_number from fields_to_plot
     if (amrex::ParallelDescriptor::NProcs() == 1){

Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp

@@ -96,6 +96,15 @@ FlushFormatCheckpoint::WriteToFile (
                      amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hzbias_fp"));
 #endif
 
+#ifdef WARPX_FERROE
+        VisMF::Write(warpx.getpolarization_fp(lev, 0),
+                     amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "px_fp"));
+        VisMF::Write(warpx.getpolarization_fp(lev, 1),
+                     amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "py_fp"));
+        VisMF::Write(warpx.getpolarization_fp(lev, 2),
+                     amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pz_fp"));
+#endif
+
         if (WarpX::fft_do_time_averaging)
         {
             VisMF::Write(warpx.getEfield_avg_fp(lev, 0),
@@ -169,6 +178,15 @@ FlushFormatCheckpoint::WriteToFile (
                          amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "Hzbias_fp"));
 #endif
 
+#ifdef WARPX_FERROE
+            VisMF::Write(warpx.getpolarization_cp(lev, 0),
+                         amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "px_cp"));
+            VisMF::Write(warpx.getpolarization_cp(lev, 1),
+                         amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "py_cp"));
+            VisMF::Write(warpx.getpolarization_cp(lev, 2),
+                         amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pz_cp"));
+#endif
+
             if (WarpX::fft_do_time_averaging)
             {
                 VisMF::Write(warpx.getEfield_avg_cp(lev, 0),

Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp

@@ -573,6 +573,11 @@ FlushFormatPlotfile::WriteAllRawFields(
         WriteRawMF( warpx.getMfield_fp(lev, 0), dm, raw_pltname, default_level_prefix, "M_xface_fp", lev, plot_raw_fields_guards);
         WriteRawMF( warpx.getMfield_fp(lev, 1), dm, raw_pltname, default_level_prefix, "M_yface_fp", lev, plot_raw_fields_guards);
         WriteRawMF( warpx.getMfield_fp(lev, 2), dm, raw_pltname, default_level_prefix, "M_zface_fp", lev, plot_raw_fields_guards);
+#endif
+#ifdef WARPX_FERROE
+        WriteRawMF( warpx.getpolarization_fp(lev, 0), dm, raw_pltname, default_level_prefix, "px_fp", lev, plot_raw_fields_guards);
+        WriteRawMF( warpx.getpolarization_fp(lev, 1), dm, raw_pltname, default_level_prefix, "py_fp", lev, plot_raw_fields_guards);
+        WriteRawMF( warpx.getpolarization_fp(lev, 2), dm, raw_pltname, default_level_prefix, "pz_fp", lev, plot_raw_fields_guards);
 #endif
         if (warpx.get_pointer_F_fp(lev))
         {

Source/Diagnostics/FullDiagnostics.cpp

@@ -16,6 +16,7 @@
 #include "Utils/WarpXAlgorithmSelection.H"
 #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
 #include "FieldSolver/London/London.H"
+#include "FieldSolver/FerroE/FerroE.H"
 #include "WarpX.H"
 
 
@@ -701,6 +702,16 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
         } else if (m_varnames[comp] == "mag_anisotropy_zface" ){
             MacroscopicProperties& macroscopic = warpx.GetMacroscopicProperties();
             m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(macroscopic.getmag_pointer_anisotropy(2), lev, m_crse_ratio);
+#endif
+#ifdef WARPX_FERROE
+        } else if ( m_varnames[comp] == "ferroelectric") {
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFerroE().m_ferroelectric_mf.get(), lev, m_crse_ratio);
+        } else if ( m_varnames[comp] == "px" ){
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_polarization_fp(lev, 0), lev, m_crse_ratio);
+        } else if ( m_varnames[comp] == "py" ){
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_polarization_fp(lev, 1), lev, m_crse_ratio);
+        } else if ( m_varnames[comp] == "pz" ){
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_polarization_fp(lev, 2), lev, m_crse_ratio);
 #endif
         } else if ( m_varnames[comp] == "superconductor") {
             m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getLondon().m_superconductor_mf.get(), lev, m_crse_ratio);

Source/Diagnostics/WarpXIO.cpp

@@ -287,6 +287,9 @@ WarpX::InitFromCheckpoint ()
             Bfield_fp[lev][i]->setVal(0.0);
 #ifdef WARPX_MAG_LLG
             Mfield_fp[lev][i]->setVal(0.0);
+#endif
+#ifdef WARPX_FERROE
+            polarization_fp[lev][i]->setVal(0.0);
 #endif
         }
 
@@ -302,6 +305,9 @@ WarpX::InitFromCheckpoint ()
                 Bfield_cp[lev][i]->setVal(0.0);
 #ifdef WARPX_MAG_LLG
                 Mfield_cp[lev][i]->setVal(0.0);
+#endif
+#ifdef WARPX_FERROE
+            polarization_cp[lev][i]->setVal(0.0);
 #endif
             }
         }
@@ -340,6 +346,16 @@ WarpX::InitFromCheckpoint ()
         VisMF::Read(*H_biasfield_fp[lev][2],
                     amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hzbias_fp"));
 #endif
+
+#ifdef WARPX_FERROE
+        VisMF::Read(*polarization_fp[lev][0],
+                    amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "px_fp"));
+        VisMF::Read(*polarization_fp[lev][1],
+                    amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "py_fp"));
+        VisMF::Read(*polarization_fp[lev][2],
+                    amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "pz_fp"));
+#endif
+
         if (WarpX::fft_do_time_averaging)
         {
             VisMF::Read(*Efield_avg_fp[lev][0],
@@ -413,6 +429,15 @@ WarpX::InitFromCheckpoint ()
             VisMF::Read(*H_biasfield_cp[lev][2],
                         amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "Hzbias_cp"));
 #endif
+
+#ifdef WARPX_FERROE
+        VisMF::Read(*polarization_cp[lev][0],
+                    amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "px_cp"));
+        VisMF::Read(*polarization_cp[lev][1],
+                    amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "py_cp"));
+        VisMF::Read(*polarization_cp[lev][2],
+                    amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "pz_cp"));
+#endif
             if (WarpX::fft_do_time_averaging)
             {
                 VisMF::Read(*Efield_avg_cp[lev][0],

Source/Initialization/WarpXInitData.cpp

@@ -750,6 +750,15 @@ WarpX::InitLevelData (int lev, Real /*time*/)
                    ::tolower);
 #endif
 
+#ifdef WARPX_FERROE
+    pp_warpx.query("P_ext_grid_init_style", P_ext_grid_s); // user-defined initial M
+    std::transform(P_ext_grid_s.begin(),
+                   P_ext_grid_s.end(),
+                   P_ext_grid_s.begin(),
+                   ::tolower);
+
+#endif
+
     // * Functions with the string "arr" in their names get an Array of
     //   values from the given entry in the table.  The array argument is
     //   resized (if necessary) to hold all the values requested.
@@ -778,6 +787,12 @@ WarpX::InitLevelData (int lev, Real /*time*/)
     if (H_bias_ext_grid_s == "constant")
         utils::parser::getArrWithParser(pp_warpx,"H_bias_external_grid", H_bias_external_grid);
 #endif
+
+#ifdef WARPX_FERROE
+    if (P_ext_grid_s == "constant")
+        utils::parser::getArrWithParser(pp_warpx, "P_external_grid", P_external_grid);
+#endif
+
     // initialize the averaged fields only if the averaged algorithm
     // is activated ('psatd.do_time_averaging=1')
     ParmParse pp_psatd("psatd");
@@ -854,6 +869,15 @@ WarpX::InitLevelData (int lev, Real /*time*/)
            }
         }
 
+#endif
+#ifdef WARPX_FERROE
+        if (P_ext_grid_s == "constant" || P_ext_grid_s == "default"){
+
+            int nghost = 1;
+            for (int icomp = 0; icomp < 3; ++icomp){
+                polarization_fp[lev][i]->setVal(P_external_grid[icomp], icomp, 2, nghost);
+            }
+        }
 #endif
    }
 
@@ -1166,6 +1190,52 @@ WarpX::InitLevelData (int lev, Real /*time*/)
 
 #endif //closes #ifdef WARPX_MAG_LLG
 
+#ifdef WARPX_FERROE
+    if (P_ext_grid_s == "parse_p_ext_grid_function") {
+#ifdef WARPX_DIM_RZ
+        amrex::Abort("P-field parser for external fields does not work with RZ");
+#endif
+        utils::parser::Store_parserString(pp_warpx, "Px_external_grid_function(x,y,z)",
+                                                    str_Px_ext_grid_function);
+        utils::parser::Store_parserString(pp_warpx, "Py_external_grid_function(x,y,z)",
+                                                    str_Py_ext_grid_function);
+        utils::parser::Store_parserString(pp_warpx, "Pz_external_grid_function(x,y,z)",
+                                                    str_Pz_ext_grid_function);
+
+        Pxfield_parser = std::make_unique<amrex::Parser>(
+                                 utils::parser::makeParser(str_Px_ext_grid_function,{"x","y","z"}));
+        Pyfield_parser = std::make_unique<amrex::Parser>(
+                                 utils::parser::makeParser(str_Py_ext_grid_function,{"x","y","z"}));
+        Pzfield_parser = std::make_unique<amrex::Parser>(
+                                 utils::parser::makeParser(str_Pz_ext_grid_function,{"x","y","z"}));
+
+       // Initialize Mfield_fp with external function directly on the faces
+       InitializeExternalFieldsOnGridUsingParser(polarization_fp[lev][0].get(),
+                                                 polarization_fp[lev][1].get(),
+                                                 polarization_fp[lev][2].get(),
+                                                 Pxfield_parser->compile<3>(),
+                                                 Pyfield_parser->compile<3>(),
+                                                 Pzfield_parser->compile<3>(),
+                                                 m_edge_lengths[lev],
+                                                 m_face_areas[lev],
+                                                 'P',
+                                                 lev, PatchType::fine);
+       if (lev > 0) {
+          InitializeExternalFieldsOnGridUsingParser(polarization_cp[lev][0].get(),
+                                                    polarization_cp[lev][1].get(),
+                                                    polarization_cp[lev][2].get(),
+                                                    Pxfield_parser->compile<3>(),
+                                                    Pyfield_parser->compile<3>(),
+                                                    Pzfield_parser->compile<3>(),
+                                                    m_edge_lengths[lev],
+                                                    m_face_areas[lev],
+                                                    'P',
+                                                    lev, PatchType::coarse);
+       }
+    }
+
+#endif //closes #ifdef WARPX_FERROE
+
     // Reading external fields from data file
     if (add_external_B_field) {
         std::string read_fields_from_path="./";

Source/WarpX.H

@@ -127,6 +127,10 @@ public:
     static amrex::Vector<amrex::Real> H_bias_external_grid;
 #endif
 
+#ifdef WARPX_FERROE
+    static amrex::Vector<amrex::Real> P_external_grid;
+#endif
+
     //! Initialization type for external magnetic field on the grid
     static std::string B_ext_grid_s;
     //! Initialization type for external electric field on the grid
@@ -145,6 +149,11 @@ public:
     static std::string M_ext_grid_s; // added for M;
     static std::string H_bias_ext_grid_s;
 #endif
+
+#ifdef WARPX_FERROE
+    static std::string P_ext_grid_s; // added for P;
+#endif
+
     static int ApplyExcitationInPML;
     //! Whether to apply the effect of an externally-defined electric field
     static bool add_external_E_field;
@@ -213,6 +222,13 @@ public:
     static std::string str_Hz_bias_ext_grid_function;
 #endif
 
+#ifdef WARPX_FERROE
+    // Parser for P_external on the grid
+    static std::string str_Px_ext_grid_function;
+    static std::string str_Py_ext_grid_function;
+    static std::string str_Pz_ext_grid_function;
+#endif
+
     //! User-defined parser to initialize x-component of the magnetic field on the grid
     std::unique_ptr<amrex::Parser> Bxfield_parser;
     //! User-defined parser to initialize y-component of the magnetic field on the grid
@@ -275,6 +291,13 @@ public:
     std::unique_ptr<amrex::Parser> Hz_biasfield_parser;
 #endif
 
+#ifdef WARPX_FERROE
+    // Parser for M_external on the grid
+    std::unique_ptr<amrex::Parser> Pxfield_parser;
+    std::unique_ptr<amrex::Parser> Pyfield_parser;
+    std::unique_ptr<amrex::Parser> Pzfield_parser;
+#endif
+
     // Algorithms
     //! Integer that corresponds to the current deposition algorithm (Esirkepov, direct, Vay)
     static short current_deposition_algo;

Source/WarpX.cpp

@@ -92,6 +92,10 @@ Vector<Real> WarpX::H_bias_external_grid(3, 0.0);
 // M could be one 9-comp vector or a vector of vectors
 #endif
 
+#ifdef WARPX_FERROE
+Vector<Real> WarpX::P_external_grid(3, 0.0);
+#endif
+
 std::string WarpX::authors = "";
 std::string WarpX::B_ext_grid_s = "default";
 std::string WarpX::E_ext_grid_s = "default";
@@ -111,6 +115,12 @@ std::string WarpX::H_ext_grid_s = "default";
 std::string WarpX::H_bias_ext_grid_s = "default";
 // "default" sets M to zero but will be overwritten by user defined input file
 #endif
+
+#ifdef WARPX_FERROE
+std::string WarpX::P_ext_grid_s = "default";
+// "default" sets P to zero but will be overwritten by user defined input file
+#endif
+
 bool WarpX::add_external_E_field = false;
 bool WarpX::add_external_B_field = false;
 
@@ -171,6 +181,13 @@ std::string WarpX::str_Hy_bias_ext_grid_function;
 std::string WarpX::str_Hz_bias_ext_grid_function;
 #endif
 
+#ifdef WARPX_FERROE
+// Parser for P_external on the grid
+std::string WarpX::str_Px_ext_grid_function;
+std::string WarpX::str_Py_ext_grid_function;
+std::string WarpX::str_Pz_ext_grid_function;
+#endif
+
 int WarpX::do_moving_window = 0;
 int WarpX::start_moving_window_step = 0;
 int WarpX::end_moving_window_step = -1;
@@ -1052,6 +1069,11 @@ WarpX::ReadParameters ()
         pp_warpx.query("mag_LLG_anisotropy_coupling",mag_LLG_anisotropy_coupling);
 #endif
 
+#ifdef WARPX_FERROE
+        // include Landau free energy contribution to effective field
+        pp_warpx.query("include_Landau", include_Landau);
+#endif
+
 #ifdef WARPX_DIM_RZ
         WARPX_ALWAYS_ASSERT_WITH_MESSAGE( isAnyBoundaryPML() == false || electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD,
             "PML are not implemented in RZ geometry with FDTD; please set a different boundary condition using boundary.field_lo and boundary.field_hi.");