diff --git a/Examples/Tests/circuits/Candice/test_metal_quibit_circuit_8_qubit_v2 b/Examples/Tests/circuits/Candice/test_metal_quibit_circuit_8_qubit_v2 new file mode 100644 index 000000000..585ae4953 --- /dev/null +++ b/Examples/Tests/circuits/Candice/test_metal_quibit_circuit_8_qubit_v2 @@ -0,0 +1,162 @@ +# See parameters in https://github.com/ECP-WarpX/artemis/pull/43 + +################################ +####### GENERAL PARAMETERS ###### +################################# +max_step = 10 + +amr.n_cell = n_cellx n_celly n_cellz +amr.max_grid_size = max_grid_sizex max_grid_sizey max_grid_sizez +amr.blocking_factor = blocking_factor +amr.refine_grid_layout = 1 # if n_MPI > n_grids, the grids will be successively divided in half until n_MPI <= n_grids + +geometry.dims = 3 +geometry.prob_lo = prob_lox prob_loy prob_loz +geometry.prob_hi = prob_hix prob_hiy prob_hiz + +amr.max_level = 0 + +# use pec instead of pml overlaying current source so you don't get a reflection +boundary.field_lo = pml pml pml +boundary.field_hi = pml pml pml + +my_constants.offset_y = -1100.e-6 + +################################# +############ NUMERICS ########### +################################# +warpx.verbose = 1 + +warpx.cfl = 0.8 + +# vacuum or macroscopic +algo.em_solver_medium = macroscopic + +# laxwendroff or backwardeuler +algo.macroscopic_sigma_method = laxwendroff + +########################################## +# sigma + +# sigma is still parsed in the PML region +macroscopic.sigma_function(x,y,z) = "sigma_0 + (sigma_si - sigma_0) * (z <= h_si)" + +# GDS file for sigma layer +macroscopic.sigma_npy_file = "/pscratch/sd/y/ytang4/Artemis_solver_veri_gds_2d/artemis/Bin/array_3d_gds_new_large_8_qubit.npy" + +# k index of where to deposit conductivity from sigma npy file +macroscopic.npy_k_index = 10 + +# enable use of PEC mask layer from sigma in valid region +algo.use_PEC_mask = 0 + +########################################## + +macroscopic.epsilon_function(x,y,z) = "eps_0 + eps_0 * (eps_r_si - 1.) * (z <= h_si)" + +macroscopic.mu_function(x,y,z) = "mu_0 + mu_0 * (mu_r_si - 1.) * (z <= h_si)" + +################################# +############ FIELDS ############# +################################# + +############### +# domain size +# n_cellx/y/z and Lx/y/z are needed to calculate dx/dy/dz +############### +my_constants.n_cellx = 3800 +my_constants.n_celly = 3800 +my_constants.n_cellz = 20 + +my_constants.max_grid_sizex = 3800 +my_constants.max_grid_sizey = 3800 +my_constants.max_grid_sizez = 20 +my_constants.blocking_factor = 2 + +my_constants.prob_lox = 0. +my_constants.prob_loy = 0. +my_constants.prob_loz = 0. + +my_constants.prob_hix = 9500.e-6 +my_constants.prob_hiy = 9500.e-6 +my_constants.prob_hiz = 200.e-6 + +my_constants.Lx = prob_hix - prob_lox +my_constants.Ly = prob_hiy - prob_loy +my_constants.Lz = prob_hiz - prob_loz + +############### +# material properties +############### +my_constants.sigma_0 = 0.0 +my_constants.sigma_si = 0. +my_constants.sigma_metal = 1.e11 + +my_constants.eps_0 = 8.8541878128e-12 +my_constants.eps_r_si = 11.7 + +my_constants.mu_0 = 1.25663706212e-06 +my_constants.mu_r_si = 1.0 +my_constants.mu_r_test = 0.9999 + +############### +# silicon and palladium cross section +############### +my_constants.h_si = 100.e-6 + +############### +# waveguide port parameters +############### +my_constants.h_port = 100.e-6 +my_constants.w_port = 200.e-6 + +############### +# derived quantities and fundamental constants - don't touch these +############### + +my_constants.pi = 3.14159265358979 + +my_constants.freq = 8.6e9 +my_constants.TP = 1./freq + +# grid spacing +my_constants.dx = Lx / n_cellx +my_constants.dy = Ly / n_celly +my_constants.dz = Lz / n_cellz + +my_constants.ddx = dx/1.e6 +my_constants.ddy = dy/1.e6 +my_constants.ddz = dz/1.e6 + +my_constants.flag_none = 0 +my_constants.flag_hs = 1 +my_constants.flag_ss = 2 + +############### +# excitation +############### + +warpx.E_excitation_on_grid_style = parse_E_excitation_grid_function + +warpx.Ex_excitation_flag_function(x,y,z) = "flag_none" +warpx.Ey_excitation_flag_function(x,y,z) = "flag_ss * ( (y >= 470.e-6 + ddy) * (y < 475.e-6 - ddy) + (y >= 485e-6 + ddy) * (y <= 490e-6 - ddy)) * (z > h_si + ddz ) * (z <= h_si + dz + ddz) * (x > 0.e-6 - ddx) * (x < 0.e-6 + ddx)" +warpx.Ez_excitation_flag_function(x,y,z) = "flag_none" + + +# This is a source on a qubit control line +warpx.Ex_excitation_grid_function(x,y,z,t) = "0." +warpx.Ey_excitation_grid_function(x,y,z,t) = "exp(-(t-3*TP)**2/(2*TP**2))*sin(2*pi*freq*t) * + ( (y >= 470.e-6 + ddy) * (y < 475.e-6 - ddy) - (y >= 485e-6 + ddy) * (y <= 490e-6 - ddy)) * (z > h_si + ddz) * (z <= h_si + dz + ddz) * (x > 0.e-6 - ddx) * (x < 0.e-6 + ddx)" +warpx.Ez_excitation_grid_function(x,y,z,t) = "0." + +############### +# diagnostics +############### + +diagnostics.diags_names = plt +############### +# full plotfiles +plt.intervals = 5 +plt.fields_to_plot = Ex Ey Ez Bx By Bz sigma mu epsilon +plt.diag_type = Full +plt.file_min_digits = 7 diff --git a/GNUmakefile b/GNUmakefile index b9f7462a0..d7fd5b6a9 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -46,3 +46,5 @@ USE_EB = FALSE WARPX_HOME := . include $(WARPX_HOME)/Source/Make.WarpX + +LIBRARIES += -lz diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 084e3f7ed..1d25078cc 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -706,10 +706,16 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri // Initialize sigma, conductivity if (macroscopic_properties->m_sigma_s == "constant") { pml_sigma_fp->setVal(macroscopic_properties->m_sigma); - } else if (macroscopic_properties->m_sigma_s == "parse_sigma_function") { + } + if (macroscopic_properties->m_sigma_s == "parse_sigma_function" || + macroscopic_properties->m_sigma_s == "parse_sigma_both") { macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_sigma_fp.get(), macroscopic_properties->m_sigma_parser->compile<3>(), lev); } + if (macroscopic_properties->m_sigma_s == "parse_sigma_npy_file" || + macroscopic_properties->m_sigma_s == "parse_sigma_both") { + macroscopic_properties->InitializeMacroMultiFabFromNumpy(pml_sigma_fp.get(), macroscopic_properties->m_sigma_npy_filename, lev, macroscopic_properties->m_npy_k_index); + } // Initialize epsilon, permittivity if (macroscopic_properties->m_epsilon_s == "constant") { @@ -875,7 +881,8 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri // Initialize sigma, conductivity if (macroscopic_properties->m_sigma_s == "constant") { pml_sigma_cp->setVal(macroscopic_properties->m_sigma); - } else if (macroscopic_properties->m_sigma_s == "parse_sigma_function") { + } else if (macroscopic_properties->m_sigma_s == "parse_sigma_function" || + macroscopic_properties->m_sigma_s == "parse_sigma_both") { macroscopic_properties->InitializeMacroMultiFabUsingParser(pml_sigma_cp.get(), macroscopic_properties->m_sigma_parser->compile<3>(), lev); } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 7820986e4..98d0fe32e 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -549,6 +549,19 @@ WarpX::OneStep_nosub (Real cur_time) ApplyExternalFieldExcitationOnGrid(ExternalFieldType::EfieldExternalPML); } + if (use_PEC_mask) { + // multiply E field by PEC mask + amrex::MultiFab& Ex = *Efield_fp[0][0].get(); + amrex::MultiFab& Ey = *Efield_fp[0][1].get(); + amrex::MultiFab& Ez = *Efield_fp[0][2].get(); + amrex::MultiFab& PECx = *PEC_fp[0][0].get(); + amrex::MultiFab& PECy = *PEC_fp[0][1].get(); + amrex::MultiFab& PECz = *PEC_fp[0][2].get(); + MultiFab::Multiply(Ex, PECx, 0, 0, 1, 0); + MultiFab::Multiply(Ey, PECy, 0, 0, 1, 0); + MultiFab::Multiply(Ez, PECz, 0, 0, 1, 0); + } + EvolveF(0.5_rt * dt[0], DtType::SecondHalf); EvolveG(0.5_rt * dt[0], DtType::SecondHalf); #ifndef WARPX_MAG_LLG diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index 966e61793..6649ccc49 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -56,7 +56,19 @@ public: amrex::ParserExecutor<3> const& macro_parser, const int lev); - /** Gpu Vector with index type of the conductivity multifab */ + void InitializeMacroMultiFabFromNumpy (amrex::MultiFab* macro_mf, + const std::string& npy_filename, + const int lev, + const int m_npy_k_index_in); + + void InitializePECFromSigma (amrex::MultiFab* sigma_mf, + amrex::MultiFab* PECx, + amrex::MultiFab* PECy, + amrex::MultiFab* PECz, + const int lev, + const int m_npy_k_index_in); + + /** Gpu Vector with index type of the conductivity multifab */ amrex::GpuArray sigma_IndexType; /** Gpu Vector with index type of the permittivity multifab */ amrex::GpuArray epsilon_IndexType; @@ -87,6 +99,11 @@ public: amrex::Real m_epsilon = PhysConst::ep0; /** Permeability, mu, of the medium */ amrex::Real m_mu = PhysConst::mu0; + /** index to insert the mask */ + int m_npy_k_index = 0; + std::string m_sigma_npy_filename; + std::string m_epsilon_npy_filename; + std::string m_mu_npy_filename; /** Parser Wrappers */ std::unique_ptr m_sigma_parser; std::unique_ptr m_epsilon_parser; @@ -123,7 +140,7 @@ public: int getmag_max_iter () {return m_mag_max_iter;} amrex::Real getmag_tol () {return m_mag_tol;} - // interpolate the magnetic properties to B locations + // interpolate the magnetic properties to B locations // magnetic properties are cell nodal // B locations are face centered // iv is an IntVect with a 1 in the face direction of interest, and 0 in the others diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index 140cba2fe..010e2f92b 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -2,6 +2,7 @@ #include "Utils/Parser/ParserUtils.H" #include "Utils/TextMsg.H" +#include "Utils/cnpy.h" #include "WarpX.H" #include @@ -41,9 +42,12 @@ MacroscopicProperties::ReadParameters () // The vacuum values are used as default for the macroscopic parameters // with a warning message to the user to indicate that no value was specified. - + // Query mask index + pp_macroscopic.query("npy_k_index", m_npy_k_index); // Query input for material conductivity, sigma. bool sigma_specified = false; + bool sigma_func_specified = false; + bool sigma_npy_specified = false; if (utils::parser::queryWithParser(pp_macroscopic, "sigma", m_sigma)) { m_sigma_s = "constant"; sigma_specified = true; @@ -51,7 +55,19 @@ MacroscopicProperties::ReadParameters () if (pp_macroscopic.query("sigma_function(x,y,z)", m_str_sigma_function) ) { m_sigma_s = "parse_sigma_function"; sigma_specified = true; + sigma_func_specified = true; + } + if (pp_macroscopic.query("sigma_npy_file", m_sigma_npy_filename) ) { + m_sigma_s = "parse_sigma_npy_file"; + sigma_specified = true; + sigma_npy_specified = true; + } + if (sigma_func_specified && sigma_npy_specified) { + // initialize both later in InitData + m_sigma_s = "parse_sigma_both"; + sigma_specified = true; } + if (!sigma_specified) { std::stringstream warnMsg; warnMsg << "Material conductivity is not specified. Using default vacuum value of " << @@ -60,7 +76,7 @@ MacroscopicProperties::ReadParameters () warnMsg.str()); } // initialization of sigma (conductivity) with parser - if (m_sigma_s == "parse_sigma_function") { + if (m_sigma_s == "parse_sigma_function" || m_sigma_s == "parse_sigma_both") { utils::parser::Store_parserString( pp_macroscopic, "sigma_function(x,y,z)", m_str_sigma_function); m_sigma_parser = std::make_unique( @@ -76,6 +92,10 @@ MacroscopicProperties::ReadParameters () m_epsilon_s = "parse_epsilon_function"; epsilon_specified = true; } + if (pp_macroscopic.query("epsilon_npy_file", m_epsilon_npy_filename) ) { + m_epsilon_s = "parse_epsilon_npy_file"; + epsilon_specified = true; + } if (!epsilon_specified) { std::stringstream warnMsg; warnMsg << "Material permittivity is not specified. Using default vacuum value of " << @@ -102,6 +122,10 @@ MacroscopicProperties::ReadParameters () m_mu_s = "parse_mu_function"; mu_specified = true; } + if (pp_macroscopic.query("mu_npy_file", m_mu_npy_filename) ) { + m_mu_s = "parse_mu_npy_file"; + mu_specified = true; + } if (!mu_specified) { std::stringstream warnMsg; warnMsg << "Material permittivity is not specified. Using default vacuum value of " << @@ -218,10 +242,32 @@ MacroscopicProperties::InitData () m_sigma_mf->setVal(m_sigma); - } else if (m_sigma_s == "parse_sigma_function") { - + } + // Step 1: Initialize sigma from parse (valid and ghost region) + if (m_sigma_s == "parse_sigma_function" || m_sigma_s == "parse_sigma_both") { InitializeMacroMultiFabUsingParser(m_sigma_mf.get(), m_sigma_parser->compile<3>(), lev); } + + // Step 2: Overwrite with numpy mask in valid region if provided + if (m_sigma_s == "parse_sigma_npy_file" || m_sigma_s == "parse_sigma_both") { + InitializeMacroMultiFabFromNumpy(m_sigma_mf.get(), m_sigma_npy_filename, lev, m_npy_k_index); + } + + if (warpx.use_PEC_mask) { + // get pointer to PEC mask in WarpX class + amrex::MultiFab * PECx = warpx.get_pointer_PEC_fp(lev,0); + amrex::MultiFab * PECy = warpx.get_pointer_PEC_fp(lev,1); + amrex::MultiFab * PECz = warpx.get_pointer_PEC_fp(lev,2); + + PECx->setVal(1.); + PECy->setVal(1.); + PECz->setVal(1.); + + InitializePECFromSigma(m_sigma_mf.get(), PECx, PECy, PECz, lev, m_npy_k_index); + // no need for sigma anymore + m_sigma_mf->setVal(0.); + } + // Initialize epsilon if (m_epsilon_s == "constant") { @@ -231,6 +277,8 @@ MacroscopicProperties::InitData () InitializeMacroMultiFabUsingParser(m_eps_mf.get(), m_epsilon_parser->compile<3>(), lev); + } else if (m_epsilon_s == "parse_epsilon_npy_file"){ + InitializeMacroMultiFabFromNumpy(m_eps_mf.get(), m_epsilon_npy_filename, lev, m_npy_k_index); } // Initialize mu @@ -242,6 +290,8 @@ MacroscopicProperties::InitData () InitializeMacroMultiFabUsingParser(m_mu_mf.get(), m_mu_parser->compile<3>(), lev); + } else if (m_mu_s == "parse_mu_npy_file"){ + InitializeMacroMultiFabFromNumpy(m_mu_mf.get(), m_mu_npy_filename, lev, m_npy_k_index); } #ifdef WARPX_MAG_LLG @@ -444,3 +494,142 @@ MacroscopicProperties::InitializeMacroMultiFabUsingParser ( } } + +void +MacroscopicProperties::InitializeMacroMultiFabFromNumpy ( + amrex::MultiFab* macro_mf, + const std::string& npy_filename, + const int lev, + const int m_npy_k_index_in) +{ + + // Load numpy array + cnpy::NpyArray arr = cnpy::npy_load(npy_filename); + amrex::Print() << "Loaded numpy array: " << npy_filename << "\n"; + + // Check expected dimensions (assume 3D) + if (arr.shape.size() != 3) { + amrex::Abort("Expected a 3D numpy array for initialization."); + } + + const size_t nx = arr.shape[0]; + const size_t ny = arr.shape[1]; + const size_t nz = arr.shape[2]; + + // Validate array size against domain + WarpX& warpx = WarpX::GetInstance(); + const amrex::Box domain_box = warpx.Geom(lev).Domain(); + const amrex::IntVect dom_size = domain_box.size(); + + // if (dom_size[0] != static_cast(nx) || + // dom_size[1] != static_cast(ny) || + // dom_size[2] != static_cast(nz)) { + // amrex::Abort("Mismatch between numpy array shape and simulation domain size!"); + // } + + if (dom_size[0] != static_cast(nx) || + dom_size[1] != static_cast(ny)) { + amrex::Abort("Mismatch between numpy array shape and simulation domain size!"); + } + + + // const double* data = arr.data(); + const double* data = arr.data(); + std::size_t total_size = nx * ny * nz; + + // Convert to amrex::Real on host + std::vector host_data(total_size); + for (std::size_t i = 0; i < total_size; ++i) { + host_data[i] = static_cast(data[i]); + } + + // Copy to device + amrex::Gpu::DeviceVector data_dev; + data_dev.resize(total_size); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, host_data.begin(), host_data.end(), data_dev.begin()); + const amrex::Real* dptr = data_dev.data(); + + amrex::IntVect iv = macro_mf->ixType().toIntVect(); + + for (MFIter mfi(*macro_mf); mfi.isValid(); ++mfi) { + // we need to loop over ghost cells since the PML algorithm creates temporary grids with valid + // and ghost regions; some may lay on the interior domain and some may lie outside the domain + const Box& bx = mfi.tilebox( iv, macro_mf->nGrowVect()); + + /* + if (mfi.index() == 0) { + const auto lo = bx.smallEnd(); + const auto hi = bx.bigEnd(); + const auto size = bx.size(); + + amrex::Print() << "fab shape: [" + << size[0] << ", " + << size[1] << ", " + << size[2] << "], index range: [" + << lo[0] << ":" << hi[0] << ", " + << lo[1] << ":" << hi[1] << ", " + << lo[2] << ":" << hi[2] << "]\n"; + } + */ + + Array4 fab = macro_mf->array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + // only overwrite macro value from the gds file if you are in the valid region + if (i >= 0 && j >= 0 && k >= 0 && i < dom_size[0] && j < dom_size[1] && k < dom_size[2]) { + if (k == m_npy_k_index_in) { + size_t idx2d = i * ny + j; // equivalent to i*ny*nz + j*nz + 0 when nz=1 + fab(i, j, k) = dptr[idx2d]; + } + } + + }); + } + amrex::Gpu::streamSynchronize(); +} + +void +MacroscopicProperties::InitializePECFromSigma (amrex::MultiFab* sigma_mf, + amrex::MultiFab* PECx, + amrex::MultiFab* PECy, + amrex::MultiFab* PECz, + const int lev, + const int m_npy_k_index_in) +{ + // PEC for Ex is on yz edges + for (MFIter mfi(*PECx); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + + Array4 sigma = sigma_mf->array(mfi); + Array4 Px = PECx->array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (k == m_npy_k_index_in) { + if (sigma(i,j,k) != 0. || sigma(i,j-1,k) != 0) { + Px(i,j,k) = 0.; + } + } + + }); + } + + // PEC for Ey is on xz edges + for (MFIter mfi(*PECy); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + + Array4 sigma = sigma_mf->array(mfi); + Array4 Py = PECy->array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (k == m_npy_k_index_in) { + if (sigma(i,j,k) != 0. || sigma(i-1,j,k) != 0) { + Py(i,j,k) = 0.; + } + } + + }); + } +} diff --git a/Source/Utils/CMakeLists.txt b/Source/Utils/CMakeLists.txt index de916cfb5..c524568b7 100644 --- a/Source/Utils/CMakeLists.txt +++ b/Source/Utils/CMakeLists.txt @@ -10,6 +10,7 @@ target_sources(WarpX WarpXUtil.cpp WarpXrocfftUtil.cpp WarpXVersion.cpp + cnpy.cpp ) add_subdirectory(Logo) diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index 83ce8543b..6fcce9557 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -9,6 +9,7 @@ CEXE_sources += IntervalsParser.cpp CEXE_sources += MPIInitHelpers.cpp CEXE_sources += RelativeCellPosition.cpp CEXE_sources += ParticleUtils.cpp +CEXE_sources += cnpy.cpp include $(WARPX_HOME)/Source/Utils/Algorithms/Make.package include $(WARPX_HOME)/Source/Utils/Logo/Make.package diff --git a/Source/Utils/cnpy.cpp b/Source/Utils/cnpy.cpp new file mode 100644 index 000000000..a2108f525 --- /dev/null +++ b/Source/Utils/cnpy.cpp @@ -0,0 +1,340 @@ +//Copyright (C) 2011 Carl Rogers +//Released under MIT License +//license available in LICENSE file, or at http://www.opensource.org/licenses/mit-license.php + +#include"cnpy.h" +#include +#include +#include +#include +#include +#include +#include +#include + +char cnpy::BigEndianTest() { + int x = 1; + return (((char *)&x)[0]) ? '<' : '>'; +} + +char cnpy::map_type(const std::type_info& t) +{ + if(t == typeid(float) ) return 'f'; + if(t == typeid(double) ) return 'f'; + if(t == typeid(long double) ) return 'f'; + + if(t == typeid(int) ) return 'i'; + if(t == typeid(char) ) return 'i'; + if(t == typeid(short) ) return 'i'; + if(t == typeid(long) ) return 'i'; + if(t == typeid(long long) ) return 'i'; + + if(t == typeid(unsigned char) ) return 'u'; + if(t == typeid(unsigned short) ) return 'u'; + if(t == typeid(unsigned long) ) return 'u'; + if(t == typeid(unsigned long long) ) return 'u'; + if(t == typeid(unsigned int) ) return 'u'; + + if(t == typeid(bool) ) return 'b'; + + if(t == typeid(std::complex) ) return 'c'; + if(t == typeid(std::complex) ) return 'c'; + if(t == typeid(std::complex) ) return 'c'; + + else return '?'; +} + +template<> std::vector& cnpy::operator+=(std::vector& lhs, const std::string rhs) { + lhs.insert(lhs.end(),rhs.begin(),rhs.end()); + return lhs; +} + +template<> std::vector& cnpy::operator+=(std::vector& lhs, const char* rhs) { + //write in little endian + size_t len = strlen(rhs); + lhs.reserve(len); + for(size_t byte = 0; byte < len; byte++) { + lhs.push_back(rhs[byte]); + } + return lhs; +} + +void cnpy::parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector& shape, bool& fortran_order) { + //std::string magic_string(buffer,6); + uint8_t major_version = *reinterpret_cast(buffer+6); + uint8_t minor_version = *reinterpret_cast(buffer+7); + uint16_t header_len = *reinterpret_cast(buffer+8); + std::string header(reinterpret_cast(buffer+9),header_len); + + size_t loc1, loc2; + + //fortran order + loc1 = header.find("fortran_order")+16; + fortran_order = (header.substr(loc1,4) == "True" ? true : false); + + //shape + loc1 = header.find("("); + loc2 = header.find(")"); + + std::regex num_regex("[0-9][0-9]*"); + std::smatch sm; + shape.clear(); + + std::string str_shape = header.substr(loc1+1,loc2-loc1-1); + while(std::regex_search(str_shape, sm, num_regex)) { + shape.push_back(std::stoi(sm[0].str())); + str_shape = sm.suffix().str(); + } + + //endian, word size, data type + //byte order code | stands for not applicable. + //not sure when this applies except for byte array + loc1 = header.find("descr")+9; + bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false); + assert(littleEndian); + + //char type = header[loc1+1]; + //assert(type == map_type(T)); + + std::string str_ws = header.substr(loc1+2); + loc2 = str_ws.find("'"); + word_size = atoi(str_ws.substr(0,loc2).c_str()); +} + +void cnpy::parse_npy_header(FILE* fp, size_t& word_size, std::vector& shape, bool& fortran_order) { + char buffer[256]; + size_t res = fread(buffer,sizeof(char),11,fp); + if(res != 11) + throw std::runtime_error("parse_npy_header: failed fread"); + std::string header = fgets(buffer,256,fp); + assert(header[header.size()-1] == '\n'); + + size_t loc1, loc2; + + //fortran order + loc1 = header.find("fortran_order"); + if (loc1 == std::string::npos) + throw std::runtime_error("parse_npy_header: failed to find header keyword: 'fortran_order'"); + loc1 += 16; + fortran_order = (header.substr(loc1,4) == "True" ? true : false); + + //shape + loc1 = header.find("("); + loc2 = header.find(")"); + if (loc1 == std::string::npos || loc2 == std::string::npos) + throw std::runtime_error("parse_npy_header: failed to find header keyword: '(' or ')'"); + + std::regex num_regex("[0-9][0-9]*"); + std::smatch sm; + shape.clear(); + + std::string str_shape = header.substr(loc1+1,loc2-loc1-1); + while(std::regex_search(str_shape, sm, num_regex)) { + shape.push_back(std::stoi(sm[0].str())); + str_shape = sm.suffix().str(); + } + + //endian, word size, data type + //byte order code | stands for not applicable. + //not sure when this applies except for byte array + loc1 = header.find("descr"); + if (loc1 == std::string::npos) + throw std::runtime_error("parse_npy_header: failed to find header keyword: 'descr'"); + loc1 += 9; + bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false); + assert(littleEndian); + + //char type = header[loc1+1]; + //assert(type == map_type(T)); + + std::string str_ws = header.substr(loc1+2); + loc2 = str_ws.find("'"); + word_size = atoi(str_ws.substr(0,loc2).c_str()); +} + +void cnpy::parse_zip_footer(FILE* fp, uint16_t& nrecs, size_t& global_header_size, size_t& global_header_offset) +{ + std::vector footer(22); + fseek(fp,-22,SEEK_END); + size_t res = fread(&footer[0],sizeof(char),22,fp); + if(res != 22) + throw std::runtime_error("parse_zip_footer: failed fread"); + + uint16_t disk_no, disk_start, nrecs_on_disk, comment_len; + disk_no = *(uint16_t*) &footer[4]; + disk_start = *(uint16_t*) &footer[6]; + nrecs_on_disk = *(uint16_t*) &footer[8]; + nrecs = *(uint16_t*) &footer[10]; + global_header_size = *(uint32_t*) &footer[12]; + global_header_offset = *(uint32_t*) &footer[16]; + comment_len = *(uint16_t*) &footer[20]; + + assert(disk_no == 0); + assert(disk_start == 0); + assert(nrecs_on_disk == nrecs); + assert(comment_len == 0); +} + +cnpy::NpyArray load_the_npy_file(FILE* fp) { + std::vector shape; + size_t word_size; + bool fortran_order; + cnpy::parse_npy_header(fp,word_size,shape,fortran_order); + + cnpy::NpyArray arr(shape, word_size, fortran_order); + size_t nread = fread(arr.data(),1,arr.num_bytes(),fp); + if(nread != arr.num_bytes()) + throw std::runtime_error("load_the_npy_file: failed fread"); + return arr; +} + +cnpy::NpyArray load_the_npz_array(FILE* fp, uint32_t compr_bytes, uint32_t uncompr_bytes) { + + std::vector buffer_compr(compr_bytes); + std::vector buffer_uncompr(uncompr_bytes); + size_t nread = fread(&buffer_compr[0],1,compr_bytes,fp); + if(nread != compr_bytes) + throw std::runtime_error("load_the_npy_file: failed fread"); + + int err; + z_stream d_stream; + + d_stream.zalloc = Z_NULL; + d_stream.zfree = Z_NULL; + d_stream.opaque = Z_NULL; + d_stream.avail_in = 0; + d_stream.next_in = Z_NULL; + err = inflateInit2(&d_stream, -MAX_WBITS); + + d_stream.avail_in = compr_bytes; + d_stream.next_in = &buffer_compr[0]; + d_stream.avail_out = uncompr_bytes; + d_stream.next_out = &buffer_uncompr[0]; + + err = inflate(&d_stream, Z_FINISH); + err = inflateEnd(&d_stream); + + std::vector shape; + size_t word_size; + bool fortran_order; + cnpy::parse_npy_header(&buffer_uncompr[0],word_size,shape,fortran_order); + + cnpy::NpyArray array(shape, word_size, fortran_order); + + size_t offset = uncompr_bytes - array.num_bytes(); + memcpy(array.data(),&buffer_uncompr[0]+offset,array.num_bytes()); + + return array; +} + +cnpy::npz_t cnpy::npz_load(std::string fname) { + FILE* fp = fopen(fname.c_str(),"rb"); + + if(!fp) { + throw std::runtime_error("npz_load: Error! Unable to open file "+fname+"!"); + } + + cnpy::npz_t arrays; + + while(1) { + std::vector local_header(30); + size_t headerres = fread(&local_header[0],sizeof(char),30,fp); + if(headerres != 30) + throw std::runtime_error("npz_load: failed fread"); + + //if we've reached the global header, stop reading + if(local_header[2] != 0x03 || local_header[3] != 0x04) break; + + //read in the variable name + uint16_t name_len = *(uint16_t*) &local_header[26]; + std::string varname(name_len,' '); + size_t vname_res = fread(&varname[0],sizeof(char),name_len,fp); + if(vname_res != name_len) + throw std::runtime_error("npz_load: failed fread"); + + //erase the lagging .npy + varname.erase(varname.end()-4,varname.end()); + + //read in the extra field + uint16_t extra_field_len = *(uint16_t*) &local_header[28]; + if(extra_field_len > 0) { + std::vector buff(extra_field_len); + size_t efield_res = fread(&buff[0],sizeof(char),extra_field_len,fp); + if(efield_res != extra_field_len) + throw std::runtime_error("npz_load: failed fread"); + } + + uint16_t compr_method = *reinterpret_cast(&local_header[0]+8); + uint32_t compr_bytes = *reinterpret_cast(&local_header[0]+18); + uint32_t uncompr_bytes = *reinterpret_cast(&local_header[0]+22); + + if(compr_method == 0) {arrays[varname] = load_the_npy_file(fp);} + else {arrays[varname] = load_the_npz_array(fp,compr_bytes,uncompr_bytes);} + } + + fclose(fp); + return arrays; +} + +cnpy::NpyArray cnpy::npz_load(std::string fname, std::string varname) { + FILE* fp = fopen(fname.c_str(),"rb"); + + if(!fp) throw std::runtime_error("npz_load: Unable to open file "+fname); + + while(1) { + std::vector local_header(30); + size_t header_res = fread(&local_header[0],sizeof(char),30,fp); + if(header_res != 30) + throw std::runtime_error("npz_load: failed fread"); + + //if we've reached the global header, stop reading + if(local_header[2] != 0x03 || local_header[3] != 0x04) break; + + //read in the variable name + uint16_t name_len = *(uint16_t*) &local_header[26]; + std::string vname(name_len,' '); + size_t vname_res = fread(&vname[0],sizeof(char),name_len,fp); + if(vname_res != name_len) + throw std::runtime_error("npz_load: failed fread"); + vname.erase(vname.end()-4,vname.end()); //erase the lagging .npy + + //read in the extra field + uint16_t extra_field_len = *(uint16_t*) &local_header[28]; + fseek(fp,extra_field_len,SEEK_CUR); //skip past the extra field + + uint16_t compr_method = *reinterpret_cast(&local_header[0]+8); + uint32_t compr_bytes = *reinterpret_cast(&local_header[0]+18); + uint32_t uncompr_bytes = *reinterpret_cast(&local_header[0]+22); + + if(vname == varname) { + NpyArray array = (compr_method == 0) ? load_the_npy_file(fp) : load_the_npz_array(fp,compr_bytes,uncompr_bytes); + fclose(fp); + return array; + } + else { + //skip past the data + uint32_t size = *(uint32_t*) &local_header[22]; + fseek(fp,size,SEEK_CUR); + } + } + + fclose(fp); + + //if we get here, we haven't found the variable in the file + throw std::runtime_error("npz_load: Variable name "+varname+" not found in "+fname); +} + +cnpy::NpyArray cnpy::npy_load(std::string fname) { + + FILE* fp = fopen(fname.c_str(), "rb"); + + if(!fp) throw std::runtime_error("npy_load: Unable to open file "+fname); + + NpyArray arr = load_the_npy_file(fp); + + fclose(fp); + return arr; +} + + + diff --git a/Source/Utils/cnpy.h b/Source/Utils/cnpy.h new file mode 100644 index 000000000..13a0f6b29 --- /dev/null +++ b/Source/Utils/cnpy.h @@ -0,0 +1,269 @@ +//Copyright (C) 2011 Carl Rogers +//Released under MIT License +//license available in LICENSE file, or at http://www.opensource.org/licenses/mit-license.php + +#ifndef LIBCNPY_H_ +#define LIBCNPY_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cnpy { + + struct NpyArray { + NpyArray(const std::vector& _shape, size_t _word_size, bool _fortran_order) : + shape(_shape), word_size(_word_size), fortran_order(_fortran_order) + { + num_vals = 1; + for(size_t i = 0;i < shape.size();i++) num_vals *= shape[i]; + data_holder = std::shared_ptr>( + new std::vector(num_vals * word_size)); + } + + NpyArray() : shape(0), word_size(0), fortran_order(0), num_vals(0) { } + + template + T* data() { + return reinterpret_cast(&(*data_holder)[0]); + } + + template + const T* data() const { + return reinterpret_cast(&(*data_holder)[0]); + } + + template + std::vector as_vec() const { + const T* p = data(); + return std::vector(p, p+num_vals); + } + + size_t num_bytes() const { + return data_holder->size(); + } + + std::shared_ptr> data_holder; + std::vector shape; + size_t word_size; + bool fortran_order; + size_t num_vals; + }; + + using npz_t = std::map; + + char BigEndianTest(); + char map_type(const std::type_info& t); + template std::vector create_npy_header(const std::vector& shape); + void parse_npy_header(FILE* fp,size_t& word_size, std::vector& shape, bool& fortran_order); + void parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector& shape, bool& fortran_order); + void parse_zip_footer(FILE* fp, uint16_t& nrecs, size_t& global_header_size, size_t& global_header_offset); + npz_t npz_load(std::string fname); + NpyArray npz_load(std::string fname, std::string varname); + NpyArray npy_load(std::string fname); + + template std::vector& operator+=(std::vector& lhs, const T rhs) { + //write in little endian + for(size_t byte = 0; byte < sizeof(T); byte++) { + char val = *((char*)&rhs+byte); + lhs.push_back(val); + } + return lhs; + } + + template<> std::vector& operator+=(std::vector& lhs, const std::string rhs); + template<> std::vector& operator+=(std::vector& lhs, const char* rhs); + + + template void npy_save(std::string fname, const T* data, const std::vector shape, std::string mode = "w") { + FILE* fp = NULL; + std::vector true_data_shape; //if appending, the shape of existing + new data + + if(mode == "a") fp = fopen(fname.c_str(),"r+b"); + + if(fp) { + //file exists. we need to append to it. read the header, modify the array size + size_t word_size; + bool fortran_order; + parse_npy_header(fp,word_size,true_data_shape,fortran_order); + assert(!fortran_order); + + if(word_size != sizeof(T)) { + std::cout<<"libnpy error: "< header = create_npy_header(true_data_shape); + size_t nels = std::accumulate(shape.begin(),shape.end(),1,std::multiplies()); + + fseek(fp,0,SEEK_SET); + fwrite(&header[0],sizeof(char),header.size(),fp); + fseek(fp,0,SEEK_END); + fwrite(data,sizeof(T),nels,fp); + fclose(fp); + } + + template void npz_save(std::string zipname, std::string fname, const T* data, const std::vector& shape, std::string mode = "w") + { + //first, append a .npy to the fname + fname += ".npy"; + + //now, on with the show + FILE* fp = NULL; + uint16_t nrecs = 0; + size_t global_header_offset = 0; + std::vector global_header; + + if(mode == "a") fp = fopen(zipname.c_str(),"r+b"); + + if(fp) { + //zip file exists. we need to add a new npy file to it. + //first read the footer. this gives us the offset and size of the global header + //then read and store the global header. + //below, we will write the the new data at the start of the global header then append the global header and footer below it + size_t global_header_size; + parse_zip_footer(fp,nrecs,global_header_size,global_header_offset); + fseek(fp,global_header_offset,SEEK_SET); + global_header.resize(global_header_size); + size_t res = fread(&global_header[0],sizeof(char),global_header_size,fp); + if(res != global_header_size){ + throw std::runtime_error("npz_save: header read error while adding to existing zip"); + } + fseek(fp,global_header_offset,SEEK_SET); + } + else { + fp = fopen(zipname.c_str(),"wb"); + } + + std::vector npy_header = create_npy_header(shape); + + size_t nels = std::accumulate(shape.begin(),shape.end(),1,std::multiplies()); + size_t nbytes = nels*sizeof(T) + npy_header.size(); + + //get the CRC of the data to be added + uint32_t crc = crc32(0L,(uint8_t*)&npy_header[0],npy_header.size()); + crc = crc32(crc,(uint8_t*)data,nels*sizeof(T)); + + //build the local header + std::vector local_header; + local_header += "PK"; //first part of sig + local_header += (uint16_t) 0x0403; //second part of sig + local_header += (uint16_t) 20; //min version to extract + local_header += (uint16_t) 0; //general purpose bit flag + local_header += (uint16_t) 0; //compression method + local_header += (uint16_t) 0; //file last mod time + local_header += (uint16_t) 0; //file last mod date + local_header += (uint32_t) crc; //crc + local_header += (uint32_t) nbytes; //compressed size + local_header += (uint32_t) nbytes; //uncompressed size + local_header += (uint16_t) fname.size(); //fname length + local_header += (uint16_t) 0; //extra field length + local_header += fname; + + //build global header + global_header += "PK"; //first part of sig + global_header += (uint16_t) 0x0201; //second part of sig + global_header += (uint16_t) 20; //version made by + global_header.insert(global_header.end(),local_header.begin()+4,local_header.begin()+30); + global_header += (uint16_t) 0; //file comment length + global_header += (uint16_t) 0; //disk number where file starts + global_header += (uint16_t) 0; //internal file attributes + global_header += (uint32_t) 0; //external file attributes + global_header += (uint32_t) global_header_offset; //relative offset of local file header, since it begins where the global header used to begin + global_header += fname; + + //build footer + std::vector footer; + footer += "PK"; //first part of sig + footer += (uint16_t) 0x0605; //second part of sig + footer += (uint16_t) 0; //number of this disk + footer += (uint16_t) 0; //disk where footer starts + footer += (uint16_t) (nrecs+1); //number of records on this disk + footer += (uint16_t) (nrecs+1); //total number of records + footer += (uint32_t) global_header.size(); //nbytes of global headers + footer += (uint32_t) (global_header_offset + nbytes + local_header.size()); //offset of start of global headers, since global header now starts after newly written array + footer += (uint16_t) 0; //zip file comment length + + //write everything + fwrite(&local_header[0],sizeof(char),local_header.size(),fp); + fwrite(&npy_header[0],sizeof(char),npy_header.size(),fp); + fwrite(data,sizeof(T),nels,fp); + fwrite(&global_header[0],sizeof(char),global_header.size(),fp); + fwrite(&footer[0],sizeof(char),footer.size(),fp); + fclose(fp); + } + + template void npy_save(std::string fname, const std::vector data, std::string mode = "w") { + std::vector shape; + shape.push_back(data.size()); + npy_save(fname, &data[0], shape, mode); + } + + template void npz_save(std::string zipname, std::string fname, const std::vector data, std::string mode = "w") { + std::vector shape; + shape.push_back(data.size()); + npz_save(zipname, fname, &data[0], shape, mode); + } + + template std::vector create_npy_header(const std::vector& shape) { + + std::vector dict; + dict += "{'descr': '"; + dict += BigEndianTest(); + dict += map_type(typeid(T)); + dict += std::to_string(sizeof(T)); + dict += "', 'fortran_order': False, 'shape': ("; + dict += std::to_string(shape[0]); + for(size_t i = 1;i < shape.size();i++) { + dict += ", "; + dict += std::to_string(shape[i]); + } + if(shape.size() == 1) dict += ","; + dict += "), }"; + //pad with spaces so that preamble+dict is modulo 16 bytes. preamble is 10 bytes. dict needs to end with \n + int remainder = 16 - (10 + dict.size()) % 16; + dict.insert(dict.end(),remainder,' '); + dict.back() = '\n'; + + std::vector header; + header += (char) 0x93; + header += "NUMPY"; + header += (char) 0x01; //major version of numpy format + header += (char) 0x00; //minor version of numpy format + header += (uint16_t) dict.size(); + header.insert(header.end(),dict.begin(),dict.end()); + + return header; + } + + +} + +#endif diff --git a/Source/WarpX.H b/Source/WarpX.H index df00c2091..346c4c0c5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -319,6 +319,8 @@ public: static int yee_coupled_solver_algo; static int lumped_inductor_algo; + static int use_PEC_mask; + //! Integer that corresponds to the order of the PSATD solution //! (whether the PSATD equations are derived from first-order or //! second-order solution) @@ -632,6 +634,7 @@ public: amrex::MultiFab * get_pointer_H_biasfield_aux (int lev, int direction) const { return H_biasfield_aux[lev][direction].get(); } #endif amrex::MultiFab * get_pointer_Efield_fp (int lev, int direction) const { return Efield_fp[lev][direction].get(); } + amrex::MultiFab * get_pointer_PEC_fp (int lev, int direction) const { return PEC_fp[lev][direction].get(); } amrex::MultiFab * get_pointer_Bfield_fp (int lev, int direction) const { return Bfield_fp[lev][direction].get(); } amrex::MultiFab * get_pointer_Bfield_sc_fp (int lev, int direction) const { return Bfield_sc_fp[lev][direction].get(); } @@ -685,6 +688,7 @@ public: const amrex::MultiFab& getcurrent_fp (int lev, int direction) {return *current_fp[lev][direction];} const amrex::MultiFab& getEfield_fp (int lev, int direction) {return *Efield_fp[lev][direction];} + const amrex::MultiFab& getPEC_fp (int lev, int direction) {return *PEC_fp[lev][direction];} const amrex::MultiFab& getBfield_fp (int lev, int direction) {return *Bfield_fp[lev][direction];} const amrex::MultiFab& getBfield_sc_fp (int lev, int direction) {return *Bfield_sc_fp[lev][direction];} #ifdef WARPX_MAG_LLG @@ -1602,6 +1606,7 @@ private: amrex::Vector, 3 > > current_fp; amrex::Vector, 3 > > current_fp_vay; amrex::Vector, 3 > > Efield_fp; + amrex::Vector, 3 > > PEC_fp; amrex::Vector, 3 > > Bfield_fp; #ifdef WARPX_MAG_LLG amrex::Vector, 3 > > Mfield_fp; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index f4fe28ed4..7dc869ab9 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -223,6 +223,7 @@ amrex::Vector WarpX::field_boundary_hi(AMREX_SPACEDIM,0); amrex::Vector WarpX::particle_boundary_lo(AMREX_SPACEDIM,ParticleBoundaryType::Absorbing); amrex::Vector WarpX::particle_boundary_hi(AMREX_SPACEDIM,ParticleBoundaryType::Absorbing); int WarpX::yee_coupled_solver_algo; +int WarpX::use_PEC_mask = 0; int WarpX::lumped_inductor_algo; bool WarpX::do_current_centering = false; @@ -387,6 +388,7 @@ WarpX::WarpX () phi_fp.resize(nlevs_max); current_fp.resize(nlevs_max); Efield_fp.resize(nlevs_max); + PEC_fp.resize(nlevs_max); Bfield_fp.resize(nlevs_max); Bfield_sc_fp.resize(nlevs_max); #ifdef WARPX_MAG_LLG @@ -1263,6 +1265,7 @@ WarpX::ReadParameters () yee_coupled_solver_algo = GetAlgorithmInteger(pp_algo, "yee_coupled_solver"); + pp_algo.query("use_PEC_mask",use_PEC_mask); lumped_inductor_algo = GetAlgorithmInteger(pp_algo, "lumped_inductor"); // Load balancing parameters @@ -1997,6 +2000,7 @@ WarpX::ClearLevel (int lev) #endif current_fp[lev][i].reset(); Efield_fp [lev][i].reset(); + PEC_fp [lev][i].reset(); Bfield_fp [lev][i].reset(); Bfield_sc_fp [lev][i].reset(); #ifdef WARPX_MAG_LLG @@ -2322,6 +2326,10 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(Efield_fp[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, tag("Efield_fp[y]")); AllocInitMultiFab(Efield_fp[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, tag("Efield_fp[z]")); + AllocInitMultiFab(PEC_fp[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, tag("PEC_fp[x]")); + AllocInitMultiFab(PEC_fp[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, tag("PEC_fp[y]")); + AllocInitMultiFab(PEC_fp[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, tag("PEC_fp[z]")); + AllocInitMultiFab(current_fp[lev][0], amrex::convert(ba, jx_nodal_flag), dm, ncomps, ngJ, tag("current_fp[x]"), 0.0_rt); AllocInitMultiFab(current_fp[lev][1], amrex::convert(ba, jy_nodal_flag), dm, ncomps, ngJ, tag("current_fp[y]"), 0.0_rt); AllocInitMultiFab(current_fp[lev][2], amrex::convert(ba, jz_nodal_flag), dm, ncomps, ngJ, tag("current_fp[z]"), 0.0_rt); diff --git a/gds_python/gds_parsing_new_modified.py b/gds_python/gds_parsing_new_modified.py new file mode 100644 index 000000000..c117a94da --- /dev/null +++ b/gds_python/gds_parsing_new_modified.py @@ -0,0 +1,75 @@ +import numpy as np +import gdspy +from skimage.draw import polygon +import matplotlib.pyplot as plt +from util import * + + + +# Load GDS +gdsii = gdspy.GdsLibrary() +gdsii.read_gds('mask7-finalNb_new_midified2.gds') +cell = gdsii.top_level()[0] + +# Get all polygons grouped by (layer, datatype) +poly_dict = cell.get_polygons(by_spec=True) + +# Extract unique layers +unique_layers = sorted(set(layer for (layer, datatype) in poly_dict.keys())) +print(f"Found {len(unique_layers)} unique layers:") +for layer in unique_layers: + print(f" - Layer {layer}") + +# Define grid size +grid_size = (4000, 4000) +array = np.zeros(grid_size, dtype=np.uint8) + +# Define which layers to include +selected_layers = {0, 1, 2, 3} # Or choose based on what's present + + +all_polys = [] +for (layer, datatype), polygons in poly_dict.items(): + if layer in selected_layers: + all_polys.extend(polygons) + +all_coords = np.vstack(all_polys) +x_min, y_min = all_coords.min(axis=0) +x_max, y_max = all_coords.max(axis=0) + + +for (layer, datatype), polygons in poly_dict.items(): + if layer not in selected_layers: + continue + for poly in polygons: + x = (poly[:, 0] - x_min) / (x_max - x_min + 1e-12) * (grid_size[1] - 1) + y = (poly[:, 1] - y_min) / (y_max - y_min + 1e-12) * (grid_size[0] - 1) + rr, cc = polygon(y, x) + rr, cc = rr.astype(int), cc.astype(int) + rr = np.clip(rr, 0, grid_size[0] - 1) + cc = np.clip(cc, 0, grid_size[1] - 1) + array[rr, cc] = 1 + +array_T = array.T +array_T_clean = remove_outer_connected_ones(array_T) +plt.figure(figsize=(6, 6)) +plt.imshow(array_T_clean, cmap='viridis', origin='lower') +plt.colorbar(label='Occupancy') +plt.title('Rasterized GDS on Layers {1,2,3}') +plt.xlabel('y') +plt.ylabel('x') +plt.tight_layout() +plt.show() + +array_T_clean_new = array_T_clean[100:3900, 100:3900] +flipped_array = 1 - array_T_clean_new + + +# Save and plot +array_3d = np.repeat(array_T_clean_new[:, :, np.newaxis], 1, axis=2) +flipped_array_3d = np.repeat(flipped_array[:, :, np.newaxis], 1, axis=2) +epsilon = 1e10 # or any small value you choose +array_3d = array_3d.astype(np.float64) * epsilon +flipped_array_3d = flipped_array_3d.astype(np.float64) * epsilon +np.save("array_3d_gds_new_large_8_qubit.npy", array_3d) +np.save("array_3d_gds_new_large_8_qubit_flipped_v2.npy", flipped_array_3d) diff --git a/gds_python/util.py b/gds_python/util.py new file mode 100644 index 000000000..abfa98910 --- /dev/null +++ b/gds_python/util.py @@ -0,0 +1,100 @@ +import numpy as np + +def pad_to_shape(array, target_shape): + """ + Pad a 2D NumPy array with zeros to the target shape, centering the original array. + + Parameters: + array (np.ndarray): 2D input array to pad. + target_shape (tuple): Desired shape (rows, cols) as (height, width). + + Returns: + np.ndarray: Zero-padded array of shape target_shape. + """ + if len(array.shape) != 2: + raise ValueError("Only 2D arrays are supported.") + + if array.shape[0] > target_shape[0] or array.shape[1] > target_shape[1]: + raise ValueError("Target shape must be greater than or equal to array shape in both dimensions.") + + pad_x = target_shape[0] - array.shape[0] + pad_y = target_shape[1] - array.shape[1] + + pad_x_before = pad_x // 2 + pad_x_after = pad_x - pad_x_before + pad_y_before = pad_y // 2 + pad_y_after = pad_y - pad_y_before + + padded = np.pad( + array, + pad_width=((pad_x_before, pad_x_after), (pad_y_before, pad_y_after)), + mode='constant', + constant_values=0 + ) + + return padded + + +import numpy as np +from scipy.ndimage import binary_fill_holes, label, generate_binary_structure + + +def remove_outer_connected_ones(mask): + """ + Remove all 1s that are connected to the image boundary. + Retains only the 'inner' objects fully surrounded by 0s. + + Parameters: + mask (np.ndarray): 2D binary array (0 and 1) + + Returns: + np.ndarray: Cleaned binary mask + """ + # Ensure binary + mask = (mask > 0).astype(np.uint8) + + # Label connected components + structure = np.ones((3, 3), dtype=np.uint8) # 8-connectivity + labeled, num_features = label(mask, structure=structure) + + # Find labels connected to the border + border_labels = set() + border_labels.update(np.unique(labeled[0, :])) # top + border_labels.update(np.unique(labeled[-1, :])) # bottom + border_labels.update(np.unique(labeled[:, 0])) # left + border_labels.update(np.unique(labeled[:, -1])) # right + + # Remove border-connected labels + cleaned_mask = np.copy(mask) + for label_id in border_labels: + cleaned_mask[labeled == label_id] = 0 + + return cleaned_mask + +import numpy as np + +def generate_2d_coordinates(shape, x_bounds, y_bounds): + """ + Generate coordinate arrays for a 2D grid. + + Parameters: + shape (tuple): Shape of the 2D array (ny, nx) + x_bounds (tuple): (x_lo, x_hi) bounds for x-axis (length = nx) + y_bounds (tuple): (y_lo, y_hi) bounds for y-axis (length = ny) + + Returns: + coords (ndarray): Array of shape (ny, nx, 2) where coords[i, j] = [x, y] + """ + nx, ny = shape + x_lo, x_hi = x_bounds + y_lo, y_hi = y_bounds + + x = np.linspace(x_lo, x_hi, nx) + y = np.linspace(y_lo, y_hi, ny) + + X, Y = np.meshgrid(x, y, indexing='xy') + coords = np.stack((X, Y), axis=-1) + + return coords + +