diff --git a/GNUmakefile b/GNUmakefile index d7fd5b6a9..b5dcabdb1 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -1,4 +1,4 @@ -AMREX_HOME ?= ../amrex +AMREX_HOME = ../../amrex DEBUG = FALSE #DEBUG = TRUE diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index f9814c84a..489c44316 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -716,6 +716,18 @@ FullDiagnostics::InitializeFieldFunctors (int lev) m_all_field_functors[lev][comp] = std::make_unique(warpx.getInductor().m_inductor_y_mf.get(), lev, m_crse_ratio); } else if ( m_varnames[comp] == "inductorz" ){ m_all_field_functors[lev][comp] = std::make_unique(warpx.getInductor().m_inductor_z_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_phi_x" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_phi_x_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_phi_y" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_phi_y_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_phi_z" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_phi_z_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_Ic_x" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_Ic_x_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_Ic_y" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_Ic_y_mf.get(), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "josephson_Ic_z" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getJosephsonJunction().m_Ic_z_mf.get(), lev, m_crse_ratio); } else { diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 65a5564f5..a160a356b 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -154,6 +154,10 @@ WarpX::Evolve (int numsteps) m_inductor->EvolveInductorJ(-0.5_rt*dt[0]); // J^(n) to J^(n-1/2) using E^(n) FillBoundaryJ(guard_cells.ng_alloc_EB); } + if (use_josephson_junction) { + m_jj->EvolveJunctionJ(-0.5_rt*dt[0]); // phi^(0) to phi^(-1/2), add JJ J^(-1/2) + FillBoundaryJ(guard_cells.ng_alloc_EB); + } is_synchronized = false; } else { if (electrostatic_solver_id == ElectrostaticSolverAlgo::None) { @@ -431,6 +435,10 @@ WarpX::OneStep_nosub (Real cur_time) m_inductor->EvolveInductorJ(dt[0]); // J^(n-1/2) to J^(n+1/2) using E^(n) FillBoundaryJ(guard_cells.ng_alloc_EB); } + if (use_josephson_junction == 1) { + m_jj->EvolveJunctionJ(dt[0]); // phi^(n-1/2) to phi^(n+1/2), add JJ J^(n+1/2) + FillBoundaryJ(guard_cells.ng_alloc_EB); + } ExecutePythonCallback("afterdeposition"); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp index 1ff8f7f60..cfa6ec973 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp @@ -133,12 +133,37 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( amrex::ignore_unused(edge_lengths); #endif + auto & warpx = WarpX::GetInstance(); + const int lev = 0; + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + + const int use_lumped_resistor = warpx.use_lumped_resistor; + const int use_lumped_capacitor = warpx.use_lumped_capacitor; + amrex::MultiFab& sigma_mf = macroscopic_properties->getsigma_mf(); amrex::MultiFab& epsilon_mf = macroscopic_properties->getepsilon_mf(); #ifndef WARPX_MAG_LLG amrex::MultiFab& mu_mf = macroscopic_properties->getmu_mf(); #endif + amrex::MultiFab* lumped_resistor_x_mf = nullptr; + amrex::MultiFab* lumped_resistor_y_mf = nullptr; + amrex::MultiFab* lumped_resistor_z_mf = nullptr; + amrex::MultiFab* lumped_capacitor_x_mf = nullptr; + amrex::MultiFab* lumped_capacitor_y_mf = nullptr; + amrex::MultiFab* lumped_capacitor_z_mf = nullptr; + + if (use_lumped_resistor) { + lumped_resistor_x_mf = macroscopic_properties->get_pointer_lumped_resistor_x(); + lumped_resistor_y_mf = macroscopic_properties->get_pointer_lumped_resistor_y(); + lumped_resistor_z_mf = macroscopic_properties->get_pointer_lumped_resistor_z(); + } + if (use_lumped_capacitor) { + lumped_capacitor_x_mf = macroscopic_properties->get_pointer_lumped_capacitor_x(); + lumped_capacitor_y_mf = macroscopic_properties->get_pointer_lumped_capacitor_y(); + lumped_capacitor_z_mf = macroscopic_properties->get_pointer_lumped_capacitor_z(); + } + // Index type required for calling ablastr::coarsen::sample::Interp to interpolate macroscopic // properties from their respective staggering to the Ex, Ey, Ez locations amrex::GpuArray const& sigma_stag = macroscopic_properties->sigma_IndexType; @@ -180,6 +205,20 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( amrex::Array4 const& mu_arr = mu_mf.array(mfi); #endif + // Lumped element arrays (default-constructed; populated only when feature enabled). + amrex::Array4 resistor_x_arr, resistor_y_arr, resistor_z_arr; + amrex::Array4 capacitor_x_arr, capacitor_y_arr, capacitor_z_arr; + if (use_lumped_resistor) { + resistor_x_arr = lumped_resistor_x_mf->array(mfi); + resistor_y_arr = lumped_resistor_y_mf->array(mfi); + resistor_z_arr = lumped_resistor_z_mf->array(mfi); + } + if (use_lumped_capacitor) { + capacitor_x_arr = lumped_capacitor_x_mf->array(mfi); + capacitor_y_arr = lumped_capacitor_y_mf->array(mfi); + capacitor_z_arr = lumped_capacitor_z_mf->array(mfi); + } + // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); int const n_coefs_x = m_stencil_coefs_x.size(); @@ -214,18 +253,25 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( // Skip field push if this cell is fully covered by embedded boundaries if (lx(i, j, k) <= 0) return; #endif - // Interpolate conductivity, sigma, to Ex position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ex_stag, macro_cr, i, j, k, scomp); - // Interpolated permittivity, epsilon, to Ex position on the grid amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ex_stag, macro_cr, i, j, k, scomp); - amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); - amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); - Ex(i, j, k) = alpha * Ex(i, j, k) - + beta * ( - T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k,0) - + T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k,0) - ) - beta * jx(i, j, k); + // fac1 = (dt/eps)*(sigma + lumped-resistor equivalent conductance) + amrex::Real const fac1 = (dt/epsilon_interp) * + ( sigma_interp + ((use_lumped_resistor && resistor_x_arr(i,j,k)!=0._rt) + ? dx[0] / (dx[1]*dx[2]*resistor_x_arr(i,j,k)) + : 0._rt) ); + // fac2 = lumped-capacitor / natural cell capacitance + amrex::Real const fac2 = (use_lumped_capacitor && capacitor_x_arr(i,j,k)!=0._rt) + ? capacitor_x_arr(i,j,k) * dx[0] / (dx[1]*dx[2]*epsilon_interp) + : 0._rt; + amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2); + amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2); + Ex(i, j, k) = alpha_c * Ex(i, j, k) + + (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k,0) + + T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k,0)) + - (dt/epsilon_interp) * beta_c * jx(i, j, k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ @@ -238,18 +284,23 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( if (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j, k)<=0 || lz(i, j-1, k)<=0) return; #endif #endif - // Interpolate conductivity, sigma, to Ey position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ey_stag, macro_cr, i, j, k, scomp); - // Interpolated permittivity, epsilon, to Ey position on the grid amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ey_stag, macro_cr, i, j, k, scomp); - amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); - amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); - Ey(i, j, k) = alpha * Ey(i, j, k) - + beta * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0) - + T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k,0) - ) - beta * jy(i, j, k); + amrex::Real const fac1 = (dt/epsilon_interp) * + ( sigma_interp + ((use_lumped_resistor && resistor_y_arr(i,j,k)!=0._rt) + ? dx[1] / (dx[0]*dx[2]*resistor_y_arr(i,j,k)) + : 0._rt) ); + amrex::Real const fac2 = (use_lumped_capacitor && capacitor_y_arr(i,j,k)!=0._rt) + ? capacitor_y_arr(i,j,k) * dx[1] / (dx[0]*dx[2]*epsilon_interp) + : 0._rt; + amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2); + amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2); + Ey(i, j, k) = alpha_c * Ey(i, j, k) + + (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0) + + T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k,0)) + - (dt/epsilon_interp) * beta_c * jy(i, j, k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ @@ -257,18 +308,23 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( // Skip field push if this cell is fully covered by embedded boundaries if (lz(i,j,k) <= 0) return; #endif - // Interpolate conductivity, sigma, to Ez position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ez_stag, macro_cr, i, j, k, scomp); - // Interpolated permittivity, epsilon, to Ez position on the grid amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ez_stag, macro_cr, i, j, k, scomp); - amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); - amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); - Ez(i, j, k) = alpha * Ez(i, j, k) - + beta * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0) - + T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k,0) - ) - beta * jz(i, j, k); + amrex::Real const fac1 = (dt/epsilon_interp) * + ( sigma_interp + ((use_lumped_resistor && resistor_z_arr(i,j,k)!=0._rt) + ? dx[2] / (dx[0]*dx[1]*resistor_z_arr(i,j,k)) + : 0._rt) ); + amrex::Real const fac2 = (use_lumped_capacitor && capacitor_z_arr(i,j,k)!=0._rt) + ? capacitor_z_arr(i,j,k) * dx[2] / (dx[0]*dx[1]*epsilon_interp) + : 0._rt; + amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2); + amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2); + Ez(i, j, k) = alpha_c * Ez(i, j, k) + + (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0) + + T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k,0)) + - (dt/epsilon_interp) * beta_c * jz(i, j, k); } ); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index 325c995e7..6fad655c5 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -47,6 +47,22 @@ public: amrex::MultiFab& getmu_mf () {return (*m_mu_mf);} amrex::MultiFab * get_pointer_mu () {return m_mu_mf.get();} + /** return MultiFabs holding the per-edge lumped resistor values [Ohm] */ + amrex::MultiFab& getlumped_resistor_x_mf () {return (*m_lumped_resistor_x_mf);} + amrex::MultiFab * get_pointer_lumped_resistor_x () {return m_lumped_resistor_x_mf.get();} + amrex::MultiFab& getlumped_resistor_y_mf () {return (*m_lumped_resistor_y_mf);} + amrex::MultiFab * get_pointer_lumped_resistor_y () {return m_lumped_resistor_y_mf.get();} + amrex::MultiFab& getlumped_resistor_z_mf () {return (*m_lumped_resistor_z_mf);} + amrex::MultiFab * get_pointer_lumped_resistor_z () {return m_lumped_resistor_z_mf.get();} + + /** return MultiFabs holding the per-edge lumped capacitor values [Farad] */ + amrex::MultiFab& getlumped_capacitor_x_mf () {return (*m_lumped_capacitor_x_mf);} + amrex::MultiFab * get_pointer_lumped_capacitor_x () {return m_lumped_capacitor_x_mf.get();} + amrex::MultiFab& getlumped_capacitor_y_mf () {return (*m_lumped_capacitor_y_mf);} + amrex::MultiFab * get_pointer_lumped_capacitor_y () {return m_lumped_capacitor_y_mf.get();} + amrex::MultiFab& getlumped_capacitor_z_mf () {return (*m_lumped_capacitor_z_mf);} + amrex::MultiFab * get_pointer_lumped_capacitor_z () {return m_lumped_capacitor_z_mf.get();} + /** Gpu Vector with index type of coarsening ratio with default value (1,1,1) */ amrex::GpuArray macro_cr_ratio; /** Initializes the Multifabs storing macroscopic properties @@ -86,6 +102,21 @@ public: /** Gpu Vector with index type of the Bz multifab */ amrex::GpuArray Bz_IndexType; + /** Gpu Vector with index type of the lumped resistor multifab */ + amrex::GpuArray lumped_resistor_x_IndexType; + amrex::GpuArray lumped_resistor_y_IndexType; + amrex::GpuArray lumped_resistor_z_IndexType; + + /** Gpu Vector with index type of the lumped capacitor multifab */ + amrex::GpuArray lumped_capacitor_x_IndexType; + amrex::GpuArray lumped_capacitor_y_IndexType; + amrex::GpuArray lumped_capacitor_z_IndexType; + + /** Gpu Vector with index type of jx/jy/jz multifab (edge-centered) */ + amrex::GpuArray jx_IndexType; + amrex::GpuArray jy_IndexType; + amrex::GpuArray jz_IndexType; + /** Stores initialization type for conductivity : constant or parser */ std::string m_sigma_s = "constant"; /** Stores initialization type for permittivity : constant or parser */ @@ -117,6 +148,14 @@ public: std::unique_ptr m_epsilon_parser; std::unique_ptr m_mu_parser; + std::unique_ptr m_lumped_resistor_x_parser; + std::unique_ptr m_lumped_resistor_y_parser; + std::unique_ptr m_lumped_resistor_z_parser; + + std::unique_ptr m_lumped_capacitor_x_parser; + std::unique_ptr m_lumped_capacitor_y_parser; + std::unique_ptr m_lumped_capacitor_z_parser; + #ifdef WARPX_MAG_LLG /** Gpu Vector with index type of the Hx multifab */ amrex::GpuArray Hx_IndexType; @@ -330,6 +369,22 @@ private: std::string m_str_epsilon_function; std::string m_str_mu_function; + /** Multifabs for lumped resistor (per-edge Ohm values) */ + std::unique_ptr m_lumped_resistor_x_mf; + std::unique_ptr m_lumped_resistor_y_mf; + std::unique_ptr m_lumped_resistor_z_mf; + std::string m_str_lumped_resistor_x_function; + std::string m_str_lumped_resistor_y_function; + std::string m_str_lumped_resistor_z_function; + + /** Multifabs for lumped capacitor (per-edge Farad values) */ + std::unique_ptr m_lumped_capacitor_x_mf; + std::unique_ptr m_lumped_capacitor_y_mf; + std::unique_ptr m_lumped_capacitor_z_mf; + std::string m_str_lumped_capacitor_x_function; + std::string m_str_lumped_capacitor_y_function; + std::string m_str_lumped_capacitor_z_function; + }; /** @@ -362,6 +417,24 @@ struct LaxWendroffAlgo { return beta; } + /** Compact form for combined dissipation (fac1) and added capacitance (fac2). + * Used when lumped resistor/capacitor terms are present. */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real alpha_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real alpha_compact = (1._rt - 0.5_rt * fac1 + fac2)/(1._rt + 0.5_rt * fac1 + fac2); + return alpha_compact; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real beta_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real beta_compact = 1._rt / (1._rt + 0.5_rt * fac1 + fac2); + return beta_compact; + } + }; /** @@ -394,6 +467,24 @@ struct BackwardEulerAlgo { return beta; } + /** Compact form for combined dissipation (fac1) and added capacitance (fac2). + * Used when lumped resistor/capacitor terms are present. */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real alpha_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real alpha_compact = (1._rt + fac2)/(1._rt + fac1 + fac2); + return alpha_compact; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real beta_compact (amrex::Real const fac1, + amrex::Real const fac2) { + using namespace amrex; + amrex::Real beta_compact = 1._rt / (1._rt + fac1 + fac2); + return beta_compact; + } + }; #endif // WARPX_MACROSCOPIC_PROPERTIES_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index 1a1774cbb..0a22fddfb 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -164,6 +164,36 @@ MacroscopicProperties::ReadParameters () utils::parser::makeParser(m_str_mu_function,{"x","y","z"})); } + // Lumped resistor R(x,y,z) [Ohm], one parser per edge direction + if (WarpX::use_lumped_resistor == 1) { + utils::parser::Store_parserString(pp_macroscopic, "lumped_resistor_x_function(x,y,z)", m_str_lumped_resistor_x_function); + m_lumped_resistor_x_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_resistor_x_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_resistor_y_function(x,y,z)", m_str_lumped_resistor_y_function); + m_lumped_resistor_y_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_resistor_y_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_resistor_z_function(x,y,z)", m_str_lumped_resistor_z_function); + m_lumped_resistor_z_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_resistor_z_function,{"x","y","z"})); + } + + // Lumped capacitor C(x,y,z) [F], one parser per edge direction + if (WarpX::use_lumped_capacitor == 1) { + utils::parser::Store_parserString(pp_macroscopic, "lumped_capacitor_x_function(x,y,z)", m_str_lumped_capacitor_x_function); + m_lumped_capacitor_x_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_capacitor_x_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_capacitor_y_function(x,y,z)", m_str_lumped_capacitor_y_function); + m_lumped_capacitor_y_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_capacitor_y_function,{"x","y","z"})); + + utils::parser::Store_parserString(pp_macroscopic, "lumped_capacitor_z_function(x,y,z)", m_str_lumped_capacitor_z_function); + m_lumped_capacitor_z_parser = std::make_unique( + utils::parser::makeParser(m_str_lumped_capacitor_z_function,{"x","y","z"})); + } + #ifdef WARPX_MAG_LLG auto &warpx = WarpX::GetInstance(); pp_macroscopic.get("mag_Ms_init_style", m_mag_Ms_s); @@ -353,6 +383,42 @@ MacroscopicProperties::InitData () } } + // Lumped elements live on the same staggering as J (edge-centered). + amrex::IntVect jx_stag = warpx.get_pointer_current_fp(lev,0)->ixType().toIntVect(); + amrex::IntVect jy_stag = warpx.get_pointer_current_fp(lev,1)->ixType().toIntVect(); + amrex::IntVect jz_stag = warpx.get_pointer_current_fp(lev,2)->ixType().toIntVect(); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + jx_IndexType[idim] = jx_stag[idim]; + jy_IndexType[idim] = jy_stag[idim]; + jz_IndexType[idim] = jz_stag[idim]; + lumped_resistor_x_IndexType[idim] = jx_stag[idim]; + lumped_resistor_y_IndexType[idim] = jy_stag[idim]; + lumped_resistor_z_IndexType[idim] = jz_stag[idim]; + lumped_capacitor_x_IndexType[idim] = jx_stag[idim]; + lumped_capacitor_y_IndexType[idim] = jy_stag[idim]; + lumped_capacitor_z_IndexType[idim] = jz_stag[idim]; + } + + if (warpx.use_lumped_resistor == 1) { + m_lumped_resistor_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, 1, ng_EB_alloc); + m_lumped_resistor_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, 1, ng_EB_alloc); + m_lumped_resistor_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, 1, ng_EB_alloc); + + InitializeMacroMultiFabUsingParser(m_lumped_resistor_x_mf.get(), m_lumped_resistor_x_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_resistor_y_mf.get(), m_lumped_resistor_y_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_resistor_z_mf.get(), m_lumped_resistor_z_parser->compile<3>(), lev); + } + + if (warpx.use_lumped_capacitor == 1) { + m_lumped_capacitor_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, 1, ng_EB_alloc); + m_lumped_capacitor_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, 1, ng_EB_alloc); + m_lumped_capacitor_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, 1, ng_EB_alloc); + + InitializeMacroMultiFabUsingParser(m_lumped_capacitor_x_mf.get(), m_lumped_capacitor_x_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_capacitor_y_mf.get(), m_lumped_capacitor_y_parser->compile<3>(), lev); + InitializeMacroMultiFabUsingParser(m_lumped_capacitor_z_mf.get(), m_lumped_capacitor_z_parser->compile<3>(), lev); + } + #ifdef WARPX_MAG_LLG // all magnetic macroparameters are stored on faces diff --git a/Source/FieldSolver/LumpedElement/CMakeLists.txt b/Source/FieldSolver/LumpedElement/CMakeLists.txt index a1f411f4a..5d4264a2e 100644 --- a/Source/FieldSolver/LumpedElement/CMakeLists.txt +++ b/Source/FieldSolver/LumpedElement/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(WarpX PRIVATE Inductor.cpp + JosephsonJunction.cpp ) diff --git a/Source/FieldSolver/LumpedElement/JosephsonJunction.H b/Source/FieldSolver/LumpedElement/JosephsonJunction.H new file mode 100644 index 000000000..55f30fe64 --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.H @@ -0,0 +1,75 @@ +/* + * This file is part of MicroEleX + * + * License: BSD-3-Clause-LBNL + */ + + +#ifndef JOSEPHSON_JUNCTION_H +#define JOSEPHSON_JUNCTION_H + +#include +#include +#include +#include +#include +#include +#include +#include + + +/** \brief Josephson junction (nonlinear inductor). + * + * Adds the supercurrent term J_J = Ic * sin(phi) / A on every edge where + * Ic(x,y,z) != 0. The phase phi obeys the second Josephson relation + * d phi / dt = (2e / hbar) * V, V = E_edge * edge_length + * advanced with leapfrog + * phi^{n+1/2} = phi^{n-1/2} + dt * (2e/hbar) * V^n. + * + * Linear resistive shunt (R) and junction capacitance (C) of an RCSJ model are + * NOT handled in this class — set them via macroscopic.lumped_resistor_* + * and macroscopic.lumped_capacitor_* on the same cell (see MacroscopicProperties). + */ +class JosephsonJunction { + +public: + JosephsonJunction (); + + void ReadParameters (); + void InitData (); + + /** Advance phi by dt using E^n, then add Ic*sin(phi)/A to the J edges. + * Pass -0.5*dt for the leapfrog initialization push, dt for the main loop. */ + void EvolveJunctionJ (amrex::Real dt); + + void InitializeJunctionMultiFabUsingParser (amrex::MultiFab *mf, + amrex::ParserExecutor<3> const& parser, + int lev); + + // Parser string (from input file) + std::string m_str_Ic_x_function; + std::string m_str_Ic_y_function; + std::string m_str_Ic_z_function; + + std::unique_ptr m_Ic_x_parser; + std::unique_ptr m_Ic_y_parser; + std::unique_ptr m_Ic_z_parser; + + /** Critical current Ic [A] per edge. Non-zero entries mark junction edges. */ + std::unique_ptr m_Ic_x_mf; + std::unique_ptr m_Ic_y_mf; + std::unique_ptr m_Ic_z_mf; + + /** Persistent phase state phi [rad], same staggering as J. */ + std::unique_ptr m_phi_x_mf; + std::unique_ptr m_phi_y_mf; + std::unique_ptr m_phi_z_mf; + + amrex::GpuArray jx_IndexType; + amrex::GpuArray jy_IndexType; + amrex::GpuArray jz_IndexType; + +}; + + +#endif diff --git a/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp b/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp new file mode 100644 index 000000000..a275bec46 --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp @@ -0,0 +1,203 @@ +#include "JosephsonJunction.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" +#include "Utils/WarpXUtil.H" +#include "Utils/WarpXConst.H" +#include "WarpX.H" +#include +#include "Utils/Parser/IntervalsParser.H" +#include "Utils/Parser/ParserUtils.H" +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +JosephsonJunction::JosephsonJunction () +{ + amrex::Print() << " JosephsonJunction (nonlinear inductor) class is constructed \n"; + ReadParameters(); +} + +void +JosephsonJunction::ReadParameters () +{ + amrex::ParmParse pp_jj("josephson"); + + utils::parser::Store_parserString(pp_jj, "Ic_x_function(x,y,z)", m_str_Ic_x_function); + m_Ic_x_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_x_function, {"x", "y", "z"})); + + utils::parser::Store_parserString(pp_jj, "Ic_y_function(x,y,z)", m_str_Ic_y_function); + m_Ic_y_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_y_function, {"x", "y", "z"})); + + utils::parser::Store_parserString(pp_jj, "Ic_z_function(x,y,z)", m_str_Ic_z_function); + m_Ic_z_parser = std::make_unique( + utils::parser::makeParser(m_str_Ic_z_function, {"x", "y", "z"})); +} + +#if( AMREX_SPACEDIM == 3) +void +JosephsonJunction::InitData() +{ + auto& warpx = WarpX::GetInstance(); + + const int lev = 0; + amrex::BoxArray ba = warpx.boxArray(lev); + amrex::DistributionMapping dmap = warpx.DistributionMap(lev); + const amrex::IntVect ng_EB_alloc = warpx.getngEB(); + + amrex::IntVect jx_stag = warpx.get_pointer_current_fp(lev,0)->ixType().toIntVect(); + amrex::IntVect jy_stag = warpx.get_pointer_current_fp(lev,1)->ixType().toIntVect(); + amrex::IntVect jz_stag = warpx.get_pointer_current_fp(lev,2)->ixType().toIntVect(); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + jx_IndexType[idim] = jx_stag[idim]; + jy_IndexType[idim] = jy_stag[idim]; + jz_IndexType[idim] = jz_stag[idim]; + } + + const int ncomps = 1; + m_Ic_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_Ic_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_Ic_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, ncomps, ng_EB_alloc); + m_phi_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_phi_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_phi_z_mf = std::make_unique(amrex::convert(ba,jz_stag), dmap, ncomps, ng_EB_alloc); + + InitializeJunctionMultiFabUsingParser(m_Ic_x_mf.get(), m_Ic_x_parser->compile<3>(), lev); + InitializeJunctionMultiFabUsingParser(m_Ic_y_mf.get(), m_Ic_y_parser->compile<3>(), lev); + InitializeJunctionMultiFabUsingParser(m_Ic_z_mf.get(), m_Ic_z_parser->compile<3>(), lev); + + m_phi_x_mf->setVal(0.0); + m_phi_y_mf->setVal(0.0); + m_phi_z_mf->setVal(0.0); +} + +void +JosephsonJunction::EvolveJunctionJ (amrex::Real dt) +{ + using namespace amrex::literals; + + amrex::Print() << " evolve Josephson junction: advance phi, add Ic*sin(phi) to J\n"; + auto & warpx = WarpX::GetInstance(); + const int lev = 0; + + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + + // 2 e / hbar = Josephson constant ~ 3.0394e15 rad/(V*s) + constexpr amrex::Real two_e_over_hbar = + 2.0_rt * PhysConst::q_e / PhysConst::hbar; + + amrex::MultiFab * jx = warpx.get_pointer_current_fp(lev, 0); + amrex::MultiFab * jy = warpx.get_pointer_current_fp(lev, 1); + amrex::MultiFab * jz = warpx.get_pointer_current_fp(lev, 2); + + amrex::MultiFab * Ex = warpx.get_pointer_Efield_fp(lev, 0); + amrex::MultiFab * Ey = warpx.get_pointer_Efield_fp(lev, 1); + amrex::MultiFab * Ez = warpx.get_pointer_Efield_fp(lev, 2); + + for (amrex::MFIter mfi(*jx, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + amrex::Array4 const& jx_arr = jx->array(mfi); + amrex::Array4 const& jy_arr = jy->array(mfi); + amrex::Array4 const& jz_arr = jz->array(mfi); + amrex::Array4 const& Ex_arr = Ex->array(mfi); + amrex::Array4 const& Ey_arr = Ey->array(mfi); + amrex::Array4 const& Ez_arr = Ez->array(mfi); + amrex::Array4 const& Ic_x_arr = m_Ic_x_mf->array(mfi); + amrex::Array4 const& Ic_y_arr = m_Ic_y_mf->array(mfi); + amrex::Array4 const& Ic_z_arr = m_Ic_z_mf->array(mfi); + amrex::Array4 const& phi_x_arr = m_phi_x_mf->array(mfi); + amrex::Array4 const& phi_y_arr = m_phi_y_mf->array(mfi); + amrex::Array4 const& phi_z_arr = m_phi_z_mf->array(mfi); + amrex::Box const& tjx = mfi.tilebox(jx->ixType().toIntVect()); + amrex::Box const& tjy = mfi.tilebox(jy->ixType().toIntVect()); + amrex::Box const& tjz = mfi.tilebox(jz->ixType().toIntVect()); + + amrex::ParallelFor(tjx, tjy, tjz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_x_arr(i,j,k) != 0.) { + const amrex::Real V = Ex_arr(i,j,k) * dx[0]; + phi_x_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[1] * dx[2]; + jx_arr(i,j,k) += Ic_x_arr(i,j,k) * std::sin(phi_x_arr(i,j,k)) / A; + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_y_arr(i,j,k) != 0.) { + const amrex::Real V = Ey_arr(i,j,k) * dx[1]; + phi_y_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[0] * dx[2]; + jy_arr(i,j,k) += Ic_y_arr(i,j,k) * std::sin(phi_y_arr(i,j,k)) / A; + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (Ic_z_arr(i,j,k) != 0.) { + const amrex::Real V = Ez_arr(i,j,k) * dx[2]; + phi_z_arr(i,j,k) += dt * two_e_over_hbar * V; + const amrex::Real A = dx[0] * dx[1]; + jz_arr(i,j,k) += Ic_z_arr(i,j,k) * std::sin(phi_z_arr(i,j,k)) / A; + } + }); + } +} + +void +JosephsonJunction::InitializeJunctionMultiFabUsingParser (amrex::MultiFab *mf, + amrex::ParserExecutor<3> const& parser, + const int lev) +{ + using namespace amrex::literals; + + WarpX& warpx = WarpX::GetInstance(); + const amrex::GpuArray dx = warpx.Geom(lev).CellSizeArray(); + const amrex::RealBox& real_box = warpx.Geom(lev).ProbDomain(); + amrex::IntVect iv = mf->ixType().toIntVect(); + for (amrex::MFIter mfi(*mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const amrex::Box& tb = mfi.tilebox(iv, mf->nGrowVect()); + amrex::Array4 const& fab = mf->array(mfi); + amrex::ParallelFor(tb, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real fac_x = (1._rt - iv[0]) * dx[0] * 0.5_rt; + amrex::Real x = i * dx[0] + real_box.lo(0) + fac_x; + amrex::Real fac_y = (1._rt - iv[1]) * dx[1] * 0.5_rt; + amrex::Real y = j * dx[1] + real_box.lo(1) + fac_y; + amrex::Real fac_z = (1._rt - iv[2]) * dx[2] * 0.5_rt; + amrex::Real z = k * dx[2] + real_box.lo(2) + fac_z; + fab(i,j,k) = parser(x,y,z); + }); + } +} + + +#else +void +JosephsonJunction::InitData() +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} + +void +JosephsonJunction::EvolveJunctionJ (amrex::Real) +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} + +void +JosephsonJunction::InitializeJunctionMultiFabUsingParser (amrex::MultiFab *, + amrex::ParserExecutor<3> const&, + const int) +{ + amrex::Abort("JosephsonJunction only works with 3D"); +} +#endif diff --git a/Source/FieldSolver/LumpedElement/Make.package b/Source/FieldSolver/LumpedElement/Make.package index a45b200ca..087c0681a 100644 --- a/Source/FieldSolver/LumpedElement/Make.package +++ b/Source/FieldSolver/LumpedElement/Make.package @@ -1,3 +1,3 @@ CEXE_sources += Inductor.cpp - +CEXE_sources += JosephsonJunction.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/LumpedElement diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 6a32c40b8..4a792717a 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -397,6 +397,16 @@ WarpX::InitData () if (WarpX::em_solver_medium==1) { m_macroscopic_properties->InitData(); } + // Lumped-element InitData must run BEFORE diagnostics InitData so that + // their multifabs exist when CellCenterFunctor captures the pointer. + if (use_lumped_inductor) { + amrex::Print() << " calling inductor (early) \n"; + m_inductor->InitData(); + } + if (use_josephson_junction) { + amrex::Print() << " calling Josephson junction (early) \n"; + m_jj->InitData(); + } InitDiagnostics(); } else @@ -407,6 +417,15 @@ WarpX::InitData () m_macroscopic_properties->InitData(); } PostRestart(); + // Same ordering constraint as the fresh-start path. + if (use_lumped_inductor) { + amrex::Print() << " calling inductor (early, restart) \n"; + m_inductor->InitData(); + } + if (use_josephson_junction) { + amrex::Print() << " calling Josephson junction (early, restart) \n"; + m_jj->InitData(); + } reduced_diags->InitData(); multi_diags->InitData(); } @@ -427,10 +446,9 @@ WarpX::InitData () m_london->InitData(); } - if (use_lumped_inductor) { - amrex::Print() << " calling inductor \n"; - m_inductor->InitData(); - } + // NOTE: Inductor and JosephsonJunction InitData were moved earlier + // (into the if/else branches above) so that diagnostics InitData sees + // their multifabs already allocated. Do not re-call here. if (ParallelDescriptor::IOProcessor()) { std::cout << "\nGrids Summary:\n"; diff --git a/Source/WarpX.H b/Source/WarpX.H index 870c7eeb5..afb335463 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -40,6 +40,7 @@ #include "Utils/WarpXAlgorithmSelection.H" #include "FieldSolver/London/London.H" #include "FieldSolver/LumpedElement/Inductor.H" +#include "FieldSolver/LumpedElement/JosephsonJunction.H" #include #include @@ -102,6 +103,7 @@ public: MacroscopicProperties& GetMacroscopicProperties () { return *m_macroscopic_properties; } London& getLondon () { return *m_london; } Inductor& getInductor () { return *m_inductor; } + JosephsonJunction& getJosephsonJunction () { return *m_jj; } MultiDiagnostics& GetMultiDiags () {return *multi_diags;} ParticleBoundaryBuffer& GetParticleBoundaryBuffer () { return *m_particle_boundary_buffer; } @@ -320,6 +322,9 @@ public: static int use_PEC_mask; static int use_lumped_inductor; + static int use_lumped_resistor; + static int use_lumped_capacitor; + static int use_josephson_junction; //! Integer that corresponds to the order of the PSATD solution //! (whether the PSATD equations are derived from first-order or @@ -1742,6 +1747,8 @@ private: std::unique_ptr m_london; // Lumped inductor std::unique_ptr m_inductor; + // Josephson junction (nonlinear inductor) + std::unique_ptr m_jj; #ifdef WARPX_MAG_LLG // time advancement scheme of M field diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 46d425e94..82f440adb 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -225,6 +225,9 @@ amrex::Vector WarpX::particle_boundary_hi(AMREX_SPACEDIM,P int WarpX::yee_coupled_solver_algo; int WarpX::use_PEC_mask = 0; int WarpX::use_lumped_inductor = 0; +int WarpX::use_lumped_resistor = 0; +int WarpX::use_lumped_capacitor = 0; +int WarpX::use_josephson_junction = 0; bool WarpX::do_current_centering = false; @@ -495,6 +498,11 @@ WarpX::WarpX () m_inductor = std::make_unique(); } + // Josephson junction (nonlinear inductor) + if (use_josephson_junction) { + m_jj = std::make_unique(); + } + // Set default values for particle and cell weights for costs update; // Default values listed here for the case AMREX_USE_GPU are determined // from single-GPU tests on Summit. @@ -1033,6 +1041,9 @@ WarpX::ReadParameters () ); } + pp_warpx.query("use_lumped_resistor", use_lumped_resistor); + pp_warpx.query("use_lumped_capacitor", use_lumped_capacitor); + #ifdef WARPX_MAG_LLG // Read the value of the time advancement scheme of M field pp_warpx.query("mag_time_scheme_order", mag_time_scheme_order); @@ -1267,6 +1278,7 @@ WarpX::ReadParameters () pp_algo.query("use_PEC_mask",use_PEC_mask); pp_algo.query("use_lumped_inductor",use_lumped_inductor); + pp_algo.query("use_josephson_junction",use_josephson_junction); // Load balancing parameters std::vector load_balance_intervals_string_vec = {"0"};