diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index 55ab056fd..a98e3c3c8 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/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){ diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 0d453c6ea..01f203eb7 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/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), diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp index 3f7bf80f8..c740a5a64 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp +++ b/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)) { diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index e233102ef..2df885d31 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/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(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(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(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(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(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(warpx.getLondon().m_superconductor_mf.get(), lev, m_crse_ratio); diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index bd6a6ad5e..6110d68d9 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/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], diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index fb3fad13a..0dc1c12a6 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/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( + utils::parser::makeParser(str_Px_ext_grid_function,{"x","y","z"})); + Pyfield_parser = std::make_unique( + utils::parser::makeParser(str_Py_ext_grid_function,{"x","y","z"})); + Pzfield_parser = std::make_unique( + 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="./"; diff --git a/Source/WarpX.H b/Source/WarpX.H index 4492cf320..2254c8166 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -127,6 +127,10 @@ public: static amrex::Vector H_bias_external_grid; #endif +#ifdef WARPX_FERROE + static amrex::Vector 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 Bxfield_parser; //! User-defined parser to initialize y-component of the magnetic field on the grid @@ -275,6 +291,13 @@ public: std::unique_ptr Hz_biasfield_parser; #endif +#ifdef WARPX_FERROE + // Parser for M_external on the grid + std::unique_ptr Pxfield_parser; + std::unique_ptr Pyfield_parser; + std::unique_ptr Pzfield_parser; +#endif + // Algorithms //! Integer that corresponds to the current deposition algorithm (Esirkepov, direct, Vay) static short current_deposition_algo; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 06ec9537d..9c0e2addc 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -92,6 +92,10 @@ Vector 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 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.");