Plasma chamber updates

Source/FieldSolver/WarpXSolveFieldsES.cpp

@@ -16,6 +16,8 @@
 void WarpX::ComputeSpaceChargeField (bool const reset_fields)
 {
     WARPX_PROFILE("WarpX::ComputeSpaceChargeField");
+    WARPX_PROFILE_REGION("WarpX::ComputeSpaceChargeField()");
+
     using ablastr::fields::Direction;
     using warpx::fields::FieldType;
 

Source/ablastr/fields/PoissonSolver.H

@@ -51,6 +51,7 @@
 #ifdef AMREX_USE_EB
 #   include <AMReX_EBFabFactory.H>
 #endif
+#include <AMReX_GMRES_MLMG.H>
 
 #include <array>
 #include <optional>
@@ -260,6 +261,31 @@ computePhi (
 
     amrex::LPInfo info;
 
+    amrex::ParmParse pp1("lpinfo");
+
+    int aaa;
+    if (pp1.query("do_agglomeration",aaa)) {
+        info.setAgglomeration(aaa); // default true
+    }
+    if (pp1.query("agg_grid_size",aaa)) {
+        info.setAgglomerationGridSize(aaa); // default -1 (architecture-dependent setting)
+    }
+    if (pp1.query("do_consolidation",aaa)) {
+        info.setConsolidation(aaa); // default true
+    }
+    if (pp1.query("con_grid_size",aaa)) {
+        info.setConsolidationGridSize(aaa); // default -1 (architecture-dependent setting)
+    }
+    if (pp1.query("con_ratio",aaa)) {
+        info.setConsolidationRatio(aaa); // default 2
+    }
+    if (pp1.query("con_strategy",aaa)) {
+        info.setConsolidationStrategy(aaa); // default 3
+    }
+    if (pp1.query("setMaxCoarseningLevel",aaa)) {
+        info.setMaxCoarseningLevel(aaa); // default 30
+    }
+
     for (int lev=0; lev<=finest_level; lev++) {
         amrex::Array<amrex::Real,AMREX_SPACEDIM> const dx_scaled
                 {AMREX_D_DECL(geom[lev].CellSize(0)/std::sqrt(1._rt-beta_solver[0]*beta_solver[0]),
@@ -399,18 +425,83 @@ computePhi (
         mlmg.setMaxIter(max_iters);
         mlmg.setAlwaysUseBNorm((max_norm_b > 0));
 
+        amrex::ParmParse pp2("mlmg");
+        int xxx;
+        amrex::Real yyy;
+
+        int solver_type = 0; // 0 = MLMG, 1 = GMRES
+        pp2.query("solver_type",solver_type);            
+
+        if (pp2.query("setPreSmooth",xxx)) {
+            mlmg.setPreSmooth(xxx); // default 2
+        }
+        if (pp2.query("setPostSmooth",xxx)) {
+            mlmg.setPostSmooth(xxx); // default 2
+        }
+        if (pp2.query("setFinalSmooth",xxx)) {
+            mlmg.setFinalSmooth(xxx); // default 8 (when smoother is used as bottom solver)
+        }
+        if (pp2.query("setBottomSmooth",xxx)) {
+            mlmg.setBottomSmooth(xxx); // default 0
+        }
+        if (pp2.query("bottomSolver",xxx)) {
+            if (xxx == 0) {
+                mlmg.setBottomSolver(amrex::BottomSolver::Default); // default 0
+            } else if (xxx == 1) {
+                mlmg.setBottomSolver(amrex::BottomSolver::smoother);
+            } else if (xxx == 2) {
+                mlmg.setBottomSolver(amrex::BottomSolver::bicgstab);
+            } else if (xxx == 3) {
+                mlmg.setBottomSolver(amrex::BottomSolver::cg);
+            } else if (xxx == 4) {
+                mlmg.setBottomSolver(amrex::BottomSolver::bicgcg);
+            } else if (xxx == 5) {
+                mlmg.setBottomSolver(amrex::BottomSolver::cgbicg);
+            } else if (xxx == 6) {
+                mlmg.setBottomSolver(amrex::BottomSolver::hypre);
+            } else if (xxx == 7) {
+                mlmg.setBottomSolver(amrex::BottomSolver::petsc);
+            }
+        }
+        if (pp2.query("setBottomTolerance",yyy)) {
+            mlmg.setBottomTolerance(yyy); // default 1.e-4
+        }
+        if (pp2.query("setBottomToleranceAbs",yyy)) {
+            mlmg.setBottomToleranceAbs(yyy); // default -1.
+        }
+
         const int ng = int(grid_type == utils::enums::GridType::Collocated); // ghost cells
         if (ng) {
             // In this case, computeE needs to use ghost nodes data. So we
             // ask MLMG to fill BC for us after it solves the problem.
             mlmg.setFinalFillBC(true);
         }
 
+        const auto mlmg_time_beg_step = static_cast<amrex::Real>(amrex::second());
+
         // Solve Poisson equation at lev
-        mlmg.solve( {phi[lev]}, {rho[lev]},
-                     relative_tolerance, absolute_tolerance );
+        if (solver_type == 0) {
+            mlmg.solve( {phi[lev]}, {rho[lev]}, relative_tolerance, absolute_tolerance );
+        } else {
+            amrex::GMRESMLMG gmsolve(mlmg);
+            gmsolve.solve( *phi[lev], *rho[lev], relative_tolerance, absolute_tolerance );
+            amrex::Vector<amrex::MultiFab> vmf;
+            vmf.emplace_back(*phi[lev], amrex::make_alias, 0, phi[lev]->nComp());
+            linop->postSolve(vmf);
+        }
 
-        const amrex::IntVect& refratio = rel_ref_ratio.value()[lev];
+#if 0
+        amrex::VisMF::Write(*rho[lev],"rho");
+        amrex::VisMF::Write(*phi[lev],"phi");
+#endif
+
+        const auto mlmg_time_end_step = static_cast<amrex::Real>(amrex::second());
+
+        amrex::Print()<< "MLMG time = "
+                      << mlmg_time_end_step-mlmg_time_beg_step
+                      << " s\n";
+
+    const amrex::IntVect& refratio = rel_ref_ratio.value()[lev];
         const int ncomp = linop->getNComp();
 
         // needed for solving the levels by levels:
@@ -428,12 +519,53 @@ computePhi (
                                         ng);
         }
 
-        // Run additional operations, such as calculation of the E field for embedded boundaries
-        if constexpr (!std::is_same_v<T_PostPhiCalculationFunctor, std::nullopt_t>) {
-            if (post_phi_calculation.has_value()) {
-                post_phi_calculation.value()(mlmg, lev);
+        if (solver_type == 0) {
+            // Run additional operations, such as calculation of the E field for embedded boundaries
+            if constexpr (!std::is_same_v<T_PostPhiCalculationFunctor, std::nullopt_t>) {
+                if (post_phi_calculation.has_value()) {
+                    post_phi_calculation.value()(mlmg, lev);
+                }
             }
         }
+
+        // create an alieas to the WarpX Efield
+        auto & warpx = WarpX::GetInstance();
+
+        amrex::Vector< amrex::Array<amrex::MultiFab *, AMREX_SPACEDIM> > e_field;
+        e_field.push_back(
+#if defined(WARPX_DIM_1D_Z)
+            amrex::Array<amrex::MultiFab*, 1>{
+                warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{2}, lev)
+            }
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+            amrex::Array<amrex::MultiFab*, 2>{
+                warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{0}, lev),
+                warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{2}, lev)
+            }
+#elif defined(WARPX_DIM_3D)
+            amrex::Array<amrex::MultiFab *, 3>{
+                warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{0}, lev),
+                warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{1}, lev),
+                warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{2}, lev)
+            }
+#endif
+        );
+
+        if (solver_type == 1) {
+            linop->compGrad(0, e_field[0], *phi[0], amrex::MLLinOp::Location::CellCenter);
+            for(int dir=0; dir<AMREX_SPACEDIM; ++dir) {
+                e_field[lev][dir]->mult(-1.);
+            }
+        }
+
+#if 0
+            amrex::VisMF::Write(*e_field[lev][0],"Ex");
+            amrex::VisMF::Write(*e_field[lev][1],"Ey");
+#if defined(WARPX_DIM_3D)
+            amrex::VisMF::Write(*e_field[lev][2],"Ez");
+#endif
+#endif
+
         rho[lev]->mult(-ablastr::constant::SI::ep0);  // Multiply rho by epsilon again
 
     } // loop over lev(els)

Tools/machines/perlmutter-nersc/perlmutter_gpu_warpx.profile.example

@@ -1,5 +1,5 @@
 # please set your project account
-export proj=""  # change me! GPU projects must end in "..._g"
+export proj="mp111_g"  # change me! GPU projects must end in "..._g"
 
 # remembers the location of this script
 export MY_PROFILE=$(cd $(dirname $BASH_SOURCE) && pwd)"/"$(basename $BASH_SOURCE)

run_plasma/CPU/inputs.2d_1024cpu

@@ -0,0 +1,130 @@
+amrex.use_gpu_aware_mpi = 1
+
+mlmg.setPreSmooth = 2
+mlmg.setPostSmooth = 2
+
+mlmg.setFinalSmooth = 8
+mlmg.setBottomSmooth = 0
+
+mlmg.bottomSolver = 0
+mlmg.setBottomTolerance = 1.e-4
+
+lpinfo.setMaxCoarseningLevel = 30
+
+# Argon
+# We want to simulate Figs. 4 (b), 5 (a) and (b), and 7 (a) in Rauf et al. 2020.
+
+my_constants.Ngas = 3.22e20  # 100 mTorr # (m^-3)
+my_constants.Tgas = 300. # (K)
+my_constants.Te = 23212. # (K) see Chen et al. 2024
+my_constants.Nplasma = 5.e16 # (m^-3) see Chen et al. 2024
+my_constants.freq = 13.56e6 # (Hz)
+my_constants.Mion = 6.63e-26 # (kg)
+my_constants.voltage = 100. # (V)
+my_constants.clight = 3.e8 # speed of light in vacuum
+my_constants.m_e = 9.11e-31 # (kg)
+my_constants.kb = 1.38e-23 # (J/K)
+
+# amr.restart = ./diags/chk01450000
+max_step = 100 # 2000000 # 5000 RF cycles
+warpx.verbose = 1
+warpx.const_dt = 1.0/(400*freq)
+warpx.do_electrostatic = labframe
+warpx.self_fields_required_precision = 1.e-7 #
+warpx.use_filter = 0
+warpx.sort_intervals = -1
+
+amr.n_cell = 4096 2048
+amr.max_grid_size_x = 128
+amr.max_grid_size_y = 64
+amr.blocking_factor = 8
+amr.max_level = 0
+
+geometry.dims = 2
+geometry.prob_lo =  -0.1035  -0.0527  # x z  # cover complete chamber (do not exploit symmetry)
+geometry.prob_hi =  0.1035   0.0527
+boundary.field_lo = pec pec
+boundary.field_hi = pec pec
+boundary.potential_hi_x = 0.
+boundary.potential_lo_z = 0.
+boundary.potential_lo_x = 0.
+boundary.potential_hi_z = 0.
+boundary.particle_lo = reflecting reflecting
+boundary.particle_hi = reflecting reflecting
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+# EB
+my_constants.te_xmax = 0.0488
+my_constants.dx_thick = 0.0032 # dielectric thickness
+my_constants.be_xmax = 0.052
+my_constants.zhi = 0.0128
+my_constants.zlo = -0.0128
+warpx.eb_implicit_function = "min(max((zlo-z),(z-zhi)),-max((x+(-be_xmax)),-(x+be_xmax)))" # "if( ((z>zhi) or (z<zlo)) and (x<be_xmax) and (x>-be_xmax) , 1,-1 )" # ??
+warpx.eb_potential(x,y,z,t) = " sin(2*pi*freq*t)*(voltage*(z>zhi)*(x<te_xmax)*(x>-te_xmax) + voltage*(z>zhi)*(x>te_xmax)*(x<be_xmax)*(be_xmax-x)/dx_thick + voltage*(z>zhi)*(x<-te_xmax)*(x>-be_xmax)*(x+be_xmax)/dx_thick + 0.*(z<zlo)*(x<be_xmax)*(x>-be_xmax)) "
+
+particles.species_names = electrons ar_ions
+
+electrons.species_type = electron
+electrons.injection_style = nuniformpercell
+electrons.initialize_self_fields = 0
+electrons.num_particles_per_cell_each_dim = 8 8
+electrons.profile = constant
+electrons.density = Nplasma
+electrons.momentum_distribution_type = maxwell_boltzmann
+electrons.theta = (kb*Te/(m_e*clight^2))
+
+ar_ions.species_type = argon
+ar_ions.charge = q_e
+ar_ions.injection_style = nuniformpercell
+ar_ions.initialize_self_fields = 0
+ar_ions.num_particles_per_cell_each_dim = 8 8
+ar_ions.profile = constant
+ar_ions.density = Nplasma
+ar_ions.momentum_distribution_type = maxwell_boltzmann
+ar_ions.theta = (kb*Tgas/(Mion*clight^2))
+ar_ions.save_particles_at_eb = 1
+
+#collisions.collision_names = coll_elec  coll_ion
+coll_ion.type = background_mcc
+coll_ion.species = ar_ions
+coll_ion.background_density = Ngas
+coll_ion.background_temperature = Tgas
+coll_ion.scattering_processes = elastic back charge_exchange
+coll_ion.elastic_cross_section = /global/cfs/cdirs/mp111/warpx-data/MCC_cross_sections/Ar/ion_scattering.dat
+coll_ion.back_cross_section = /global/cfs/cdirs/mp111/warpx-data/MCC_cross_sections/Ar/ion_back_scatter.dat
+coll_ion.charge_exchange_cross_section = /global/cfs/cdirs/mp111/warpx-data/MCC_cross_sections/Ar/charge_exchange.dat
+
+coll_elec.type = background_mcc
+coll_elec.species = electrons
+coll_elec.background_density = Ngas
+coll_elec.background_temperature = Tgas
+coll_elec.scattering_processes = elastic excitation1 ionization
+coll_elec.elastic_cross_section = /global/cfs/cdirs/mp111/warpx-data/MCC_cross_sections/Ar/electron_scattering.dat
+coll_elec.excitation1_energy = 11.5
+coll_elec.excitation1_cross_section = /global/cfs/cdirs/mp111/warpx-data/MCC_cross_sections/Ar/excitation_1.dat
+coll_elec.ionization_energy = 15.7596112
+coll_elec.ionization_cross_section = /global/cfs/cdirs/mp111/warpx-data/MCC_cross_sections/Ar/ionization.dat
+coll_elec.ionization_species = ar_ions
+
+#diagnostics.diags_names = plt chk #plt_eb
+#plt.diag_type = Full
+#plt.intervals = 1
+#plt.fields_to_plot = phi
+#plt.file_min_digits = 8
+
+#plt_eb.diag_type = BoundaryScraping
+#plt_eb.format = openpmd
+#plt_eb.fields_to_plot = phi
+#plt_eb.particle_field_species = ar_ions
+#plt_eb.intervals = 190000:200000:1000
+
+#chk.diag_type = Full
+#chk.format = checkpoint
+#chk.intervals = 1000
+#chk.file_min_digits = 8
+
+#warpx.reduced_diags_names = partnum
+#partnum.type = ParticleNumber
+#partnum.intervals = 1