diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 249f204a9..97a681a96 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp @@ -115,8 +115,10 @@ FlushFormatCheckpoint::WriteToFile ( if (warpx.getis_synchronized() || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || warpx.use_lumped_inductor == 1) { + || warpx.use_lumped_inductor == 1 + || warpx.use_josephson_junction == 1) { // Need to save j if synchronized because after restart we need j to evolve E by dt/2. + // London / lumped inductor / Josephson junction all treat j as persistent state. VisMF::Write(warpx.getcurrent_fp(lev, 0), amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "jx_fp")); VisMF::Write(warpx.getcurrent_fp(lev, 1), @@ -180,7 +182,8 @@ FlushFormatCheckpoint::WriteToFile ( if (warpx.getis_synchronized() || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || warpx.use_lumped_inductor == 1) { + || warpx.use_lumped_inductor == 1 + || warpx.use_josephson_junction == 1) { // 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")); @@ -205,6 +208,12 @@ FlushFormatCheckpoint::WriteToFile ( } } + // Josephson junction state (phi, Js). phi is the junction's only real state and + // is NOT reconstructible from the fields, so without this a restart resets it to 0. + if (warpx.use_josephson_junction == 1) { + warpx.getJosephsonJunction().WriteCheckpoint(checkpointname); + } + CheckpointParticles(checkpointname, particle_diags); WriteDMaps(checkpointname, nlev); 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/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index ed617acaf..025ac5393 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -359,7 +359,8 @@ WarpX::InitFromCheckpoint () if (is_synchronized || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || use_lumped_inductor == 1) + || use_lumped_inductor == 1 + || use_josephson_junction == 1) { VisMF::Read(*current_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_fp")); @@ -426,7 +427,8 @@ WarpX::InitFromCheckpoint () if (is_synchronized || WarpX::yee_coupled_solver_algo == CoupledYeeSolver::MaxwellLondon - || use_lumped_inductor == 1) { + || use_lumped_inductor == 1 + || use_josephson_junction == 1) { VisMF::Read(*current_cp[lev][0], amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "jx_cp")); VisMF::Read(*current_cp[lev][1], diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 6303c0f71..348ae238f 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/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..61fa6a8ea --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.H @@ -0,0 +1,97 @@ +/* + * 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). + * + * Sets 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. + * Because the current MultiFab is a persistent state in this solver (it is not + * re-zeroed each step), the junction updates J by the *difference* between the + * new and previously-added supercurrent so that J holds the instantaneous + * Ic*sin(phi)/A, not a running sum. This also composes with any London / + * inductor current sharing the same edge. + * + * 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); + + /** Persist the junction state (phi, Js) across a restart. + * Without these, a restart silently resets phi to 0 while the EM fields are + * restored from the checkpoint, leaving the junction phase inconsistent with + * the fields (the JJ loses its history and the post-restart dynamics are wrong). + * Call WriteCheckpoint when the checkpoint is flushed, and ReadCheckpoint AFTER + * InitData() on restart (InitData allocates and zeroes these MultiFabs). */ + void WriteCheckpoint (const std::string& dir) const; + void ReadCheckpoint (const std::string& dir); + + // 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; + + /** Supercurrent density Ic*sin(phi)/A [A/m^2] added to J on the previous + * push. J is a persistent state in this solver (never re-zeroed), so we + * subtract this on the next push and add the new value — J then carries the + * instantaneous supercurrent rather than a running sum over time steps. */ + std::unique_ptr m_Js_x_mf; + std::unique_ptr m_Js_y_mf; + std::unique_ptr m_Js_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..7cc395f6f --- /dev/null +++ b/Source/FieldSolver/LumpedElement/JosephsonJunction.cpp @@ -0,0 +1,267 @@ +#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 + +#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); + m_Js_x_mf = std::make_unique(amrex::convert(ba,jx_stag), dmap, ncomps, ng_EB_alloc); + m_Js_y_mf = std::make_unique(amrex::convert(ba,jy_stag), dmap, ncomps, ng_EB_alloc); + m_Js_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); + m_Js_x_mf->setVal(0.0); + m_Js_y_mf->setVal(0.0); + m_Js_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::Array4 const& Js_x_arr = m_Js_x_mf->array(mfi); + amrex::Array4 const& Js_y_arr = m_Js_y_mf->array(mfi); + amrex::Array4 const& Js_z_arr = m_Js_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]; + const amrex::Real Js_new = Ic_x_arr(i,j,k) * std::sin(phi_x_arr(i,j,k)) / A; + jx_arr(i,j,k) += Js_new - Js_x_arr(i,j,k); + Js_x_arr(i,j,k) = Js_new; + } + }, + [=] 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]; + const amrex::Real Js_new = Ic_y_arr(i,j,k) * std::sin(phi_y_arr(i,j,k)) / A; + jy_arr(i,j,k) += Js_new - Js_y_arr(i,j,k); + Js_y_arr(i,j,k) = Js_new; + } + }, + [=] 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]; + const amrex::Real Js_new = Ic_z_arr(i,j,k) * std::sin(phi_z_arr(i,j,k)) / A; + jz_arr(i,j,k) += Js_new - Js_z_arr(i,j,k); + Js_z_arr(i,j,k) = Js_new; + } + }); + } +} + +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 + + +// --------------------------------------------------------------------------- +// Checkpoint the persistent junction state (phi, Js). +// phi is the junction's only real state; Js is the supercurrent this class +// last added to J (needed by the delta update in EvolveJunctionJ). Both are +// dimension-agnostic MultiFab I/O, so they live outside the 3D-only block. +// --------------------------------------------------------------------------- +void +JosephsonJunction::WriteCheckpoint (const std::string& dir) const +{ + if (!m_phi_x_mf) { return; } // junction never initialized + const int lev = 0; + const std::string lp{"Level_"}; + amrex::VisMF::Write(*m_phi_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_x")); + amrex::VisMF::Write(*m_phi_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_y")); + amrex::VisMF::Write(*m_phi_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_z")); + amrex::VisMF::Write(*m_Js_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_x")); + amrex::VisMF::Write(*m_Js_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_y")); + amrex::VisMF::Write(*m_Js_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_z")); +} + +void +JosephsonJunction::ReadCheckpoint (const std::string& dir) +{ + if (!m_phi_x_mf) { return; } + const int lev = 0; + const std::string lp{"Level_"}; + // Checkpoints written before the junction state was saved have no jj_* files. + // Warn and keep the freshly-zeroed phi rather than aborting -- but the restart + // is then NOT physically continuous (phi loses its history). + if (!amrex::FileExists(amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_x") + "_H")) { + amrex::Print() << "WARNING: checkpoint '" << dir << "' contains no Josephson junction " + << "state (jj_phi_*). phi restarts from 0 while the fields are restored, " + << "so the junction is inconsistent with the fields -- rerun from scratch " + << "or use a checkpoint written with JJ state.\n"; + return; + } + amrex::VisMF::Read(*m_phi_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_x")); + amrex::VisMF::Read(*m_phi_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_y")); + amrex::VisMF::Read(*m_phi_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_phi_z")); + amrex::VisMF::Read(*m_Js_x_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_x")); + amrex::VisMF::Read(*m_Js_y_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_y")); + amrex::VisMF::Read(*m_Js_z_mf, amrex::MultiFabFileFullPrefix(lev, dir, lp, "jj_Js_z")); + amrex::Print() << " Josephson junction state (phi, Js) restored from checkpoint\n"; +} 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 db63185f9..18ef97e83 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -424,6 +424,23 @@ WarpX::InitData () m_london->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 \n"; + m_inductor->InitData(); + } + if (use_josephson_junction) { + amrex::Print() << " calling Josephson junction \n"; + m_jj->InitData(); + if (!restart_chkfile.empty()) { + // InitData() allocates and ZEROES phi/Js. Restore them from the checkpoint + // so the junction continues its phase; otherwise phi restarts at 0 while the + // fields are restored, which silently corrupts the post-restart dynamics. + m_jj->ReadCheckpoint(restart_chkfile); + } + } + if (restart_chkfile.empty()) { InitDiagnostics(); } else { @@ -431,11 +448,6 @@ WarpX::InitData () multi_diags->InitData(); } - if (use_lumped_inductor) { - amrex::Print() << " calling inductor \n"; - m_inductor->InitData(); - } - if (ParallelDescriptor::IOProcessor()) { std::cout << "\nGrids Summary:\n"; printGridSummary(std::cout, 0, finestLevel()); diff --git a/Source/WarpX.H b/Source/WarpX.H index c5aa50c5e..f89a0e8b9 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; } @@ -336,6 +338,7 @@ public: 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 @@ -1758,6 +1761,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 c49e65a9f..d3abb397b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -233,6 +233,7 @@ int WarpX::extend_PEC_mask_into_pml = 1; 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; @@ -503,6 +504,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. @@ -1305,6 +1311,7 @@ WarpX::ReadParameters () pp_algo.query("use_PEC_mask",use_PEC_mask); pp_algo.query("extend_PEC_mask_into_pml",extend_PEC_mask_into_pml); 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"};