Ferrorelectric support in ARTEMIS

GNUmakefile

@@ -41,7 +41,8 @@ WarpxBinDir = Bin
 USE_PSATD = FALSE
 USE_PSATD_PICSAR = FALSE
 USE_RZ = FALSE
-USE_LLG = TRUE
+USE_LLG = FALSE
+USE_FERROE = TRUE
 
 USE_EB = FALSE
 

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),
@@ -123,6 +132,16 @@ FlushFormatCheckpoint::WriteToFile (
                          amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jz_fp"));
         }
 
+        if (warpx.getis_synchronized() || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellFerroE) {
+            // Need to save j if synchronized because after restart we need j to evolve E by dt/2.
+            VisMF::Write(warpx.getcurrent_fp(lev, 0),
+                         amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jx_fp"));
+            VisMF::Write(warpx.getcurrent_fp(lev, 1),
+                         amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jy_fp"));
+            VisMF::Write(warpx.getcurrent_fp(lev, 2),
+                         amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jz_fp"));
+        }
+
         if (lev > 0)
         {
             VisMF::Write(warpx.getEfield_cp(lev, 0),
@@ -159,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),
@@ -185,6 +213,16 @@ FlushFormatCheckpoint::WriteToFile (
                 VisMF::Write(warpx.getcurrent_cp(lev, 2),
                              amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jz_cp"));
             }
+
+            if (warpx.getis_synchronized() || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellFerroE) {
+                // Need to save j if synchronized because after restart we need j to evolve E by dt/2.
+                VisMF::Write(warpx.getcurrent_cp(lev, 0),
+                             amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jx_cp"));
+                VisMF::Write(warpx.getcurrent_cp(lev, 1),
+                             amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jy_cp"));
+                VisMF::Write(warpx.getcurrent_cp(lev, 2),
+                             amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jz_cp"));
+            }
         }
 
         if (warpx.DoPML()) {

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,22 @@ 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, true, 1, 0);
+        } 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, true, 1, 0);
+        } 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, true, 1, 0);
+        } else if ( m_varnames[comp] == "dpxdt" ){
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_polarization_fp(lev, 0), lev, m_crse_ratio, true, 1, 1);
+        } else if ( m_varnames[comp] == "dpydt" ){
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_polarization_fp(lev, 1), lev, m_crse_ratio, true, 1, 1);
+        } else if ( m_varnames[comp] == "dpzdt" ){
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_polarization_fp(lev, 2), lev, m_crse_ratio, true, 1, 1);
 #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],
@@ -366,6 +382,15 @@ WarpX::InitFromCheckpoint ()
                         amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jz_fp"));
         }
 
+        if (is_synchronized || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellFerroE) {
+            VisMF::Read(*current_fp[lev][0],
+                        amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_fp"));
+            VisMF::Read(*current_fp[lev][1],
+                        amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jy_fp"));
+            VisMF::Read(*current_fp[lev][2],
+                        amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jz_fp"));
+        }
+
         if (lev > 0)
         {
             VisMF::Read(*Efield_cp[lev][0],
@@ -404,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],
@@ -429,6 +463,15 @@ WarpX::InitFromCheckpoint ()
                 VisMF::Read(*current_cp[lev][2],
                             amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jz_cp"));
             }
+
+            if (is_synchronized || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellFerroE) {
+                VisMF::Read(*current_cp[lev][0],
+                            amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_cp"));
+                VisMF::Read(*current_cp[lev][1],
+                            amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jy_cp"));
+                VisMF::Read(*current_cp[lev][2],
+                            amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jz_cp"));
+            }
         }
     }
 

Source/Evolve/WarpXEvolve.cpp

@@ -32,6 +32,7 @@
 #include "Utils/WarpXProfilerWrapper.H"
 #include "Utils/WarpXUtil.H"
 #include "FieldSolver/London/London.H"
+#include "FieldSolver/FerroE/FerroE.H"
 
 #include <ablastr/utils/SignalHandling.H>
 #include <ablastr/warn_manager/WarnManager.H>
@@ -150,6 +151,14 @@ WarpX::Evolve (int numsteps)
                 m_london->EvolveLondonJ(-0.5_rt*dt[0]); // J^(n) to J^(n-1/2) using E^(n)
                 FillBoundaryJ(guard_cells.ng_alloc_EB);
             }
+#ifdef WARPX_FERROE
+            if (WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellFerroE) {
+                m_ferroe->EvolveP(-0.5_rt*dt[0]); // P^(n) to P^(n-1/2) using E^(n)
+                FillBoundaryP(guard_cells.ng_alloc_EB);
+                m_ferroe->EvolveFerroEJ(-0.5_rt*dt[0]); // J^(n) to J^(n-1/2) using E^(n)
+                FillBoundaryJ(guard_cells.ng_alloc_EB);
+            }
+#endif
             is_synchronized = false;
         } else {
             if (electrostatic_solver_id == ElectrostaticSolverAlgo::None) {
@@ -203,6 +212,9 @@ WarpX::Evolve (int numsteps)
             if (WarpX::yee_coupled_solver_algo != CoupledYeeSolver::MaxwellLondon) {
                 PushParticlesandDepose(cur_time, skip_deposition);
             }
+            if (WarpX::yee_coupled_solver_algo != CoupledYeeSolver::MaxwellFerroE) {
+                PushParticlesandDepose(cur_time, skip_deposition);
+            }
         }
         // Electromagnetic case: multi-J algorithm
         else if (do_multi_J)
@@ -429,6 +441,9 @@ WarpX::OneStep_nosub (Real cur_time)
     if (WarpX::yee_coupled_solver_algo != CoupledYeeSolver::MaxwellLondon) {
         PushParticlesandDepose(cur_time);
     }
+    if (WarpX::yee_coupled_solver_algo != CoupledYeeSolver::MaxwellFerroE) {
+        PushParticlesandDepose(cur_time);
+    }
 #ifndef WARPX_MAG_LLG
     if (WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon) {
         amrex::Print() << " in evolve london j\n";
@@ -438,6 +453,16 @@ WarpX::OneStep_nosub (Real cur_time)
         // fill boundary here
     }
 #endif
+#ifdef WARPX_FERROE
+    if (WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellFerroE) {
+        amrex::Print() << " in evolve ferroelectric j\n";
+        m_ferroe->EvolveP(dt[0]); // P^(n-1/2) to P^(n+1/2) using E^(n)
+        FillBoundaryP(guard_cells.ng_alloc_EB);
+        m_ferroe->EvolveFerroEJ(dt[0]); // J^(n-1/2) to J^(n+1/2) using E^(n)
+        FillBoundaryJ(guard_cells.ng_alloc_EB);
+        // fill boundary here
+    }
+#endif
 
     ExecutePythonCallback("afterdeposition");
 

Source/FieldSolver/FerroE/CMakeLists.txt

@@ -0,0 +1,4 @@
+target_sources(WarpX
+  PRIVATE
+  FerroE.cpp
+)