Pulsar remove managed memory allocation

Examples/Physics_applications/pulsar/inputs.corotating.3d.PM

@@ -13,16 +13,18 @@
 #################################
 ####### GENERAL PARAMETERS ######
 #################################
-max_step = 5
+max_step = 1
 amr.n_cell = 128 128 128
 amr.max_grid_size = 128
 amr.blocking_factor = 128
 amr.max_level = 0
 geometry.coord_sys   = 0                  # 0: Cartesian
-geometry.is_periodic = 0 0 0     # Is periodic?
-geometry.prob_lo     = 0.0    0.0    0.0
+geometry.prob_lo     = 0 0 0
 geometry.prob_hi     = 180000 180000 180000
 
+boundary.field_lo = pml pml pml
+boundary.field_hi = pml pml pml
+
 #################################
 ############ NUMERICS ###########
 #################################
@@ -31,26 +33,25 @@ warpx.verbose = 1
 warpx.do_dive_cleaning = 0
 warpx.use_filter = 1
 warpx.cfl = .95
-warpx.do_pml = 1
 my_constants.pi    = 3.141592653589793
 my_constants.dens  = 5.544e6
 my_constants.xc  = 90000
 my_constants.yc  = 90000
 my_constants.zc  = 90000
-my_constants.r_star  = 12000
+my_constants.r_star  = 12032
 my_constants.omega = 6245.676
 my_constants.c = 299792458.
 my_constants.B_star = 8.0323e-6                  # magnetic field of NS (T)
-my_constants.dR = 1400
+my_constants.dR = 1406.25
 my_constants.to = 2.e-4
 algo.particle_shape = 3
 
 
 #################################
 ############ PLASMA #############
 #################################
-particles.nspecies = 2
-particles.species_names = plasma_e plasma_p
+#particles.nspecies = 2
+#particles.species_names = plasma_e plasma_p
 
 plasma_e.charge = -q_e
 plasma_e.mass = m_e
@@ -102,11 +103,11 @@ pulsar.ramp_omega_time = 0.0             # time over which to ramp up NS angular
 pulsar.center_star = 90000 90000 90000
 pulsar.R_star = 12032                  # radius of NS (m)
 pulsar.B_star = 8.0323e-6                  # magnetic field of NS (T)
-pulsar.dR = 1400                       # thickness on the surface of the pulsar
+pulsar.dR = 1406.25                       # thickness on the surface of the pulsar
                                        # consistency requires dR = my_constants.dR
 pulsar.verbose = 0                     # [0/1]: turn on verbosity for debugging print statements
 pulsar.EB_external = 0                 # [0/1]: to apply external E and B
-pulsar.E_external_monopole = 0                  # [0/1]
+pulsar.E_external_monopole = 1                  # [0/1]
 pulsar.EB_corotating_maxradius = 12032 # The radius where E-field changes from corotating
                                        # inside the star to quadrapole outside.
                                        # Default is R_star. (r<=cor_radius) i.e. the includes
@@ -135,11 +136,18 @@ pulsar.turnoff_plasmaEB_gather = 0    # (0/1) 0 is default. always gather
 pulsar.max_nogather_radius = 0.       # radius within which self-consistent field gather
 pulsar.init_dipoleBfield = 1          # default is 1
 pulsar.init_corotatingEfield = 1      # default
-pulsar.corotatingE_maxradius = 12032  # default is Rstar
+pulsar.corotatingE_maxradius = 16250.75  # default is Rstar
 pulsar.enforceDipoleB_maxradius = 12032 # default is Rstar
 pulsar.enforceCorotatingE = 1        # default 1
 pulsar.enforceDipoleB = 1        # default 1
 pulsar.singleParticleTest = 0    # default 0, if 1, then a particle pair will be introduced.
                                  # Single particle position/momentum will need to be specified
-pulsar.DampBDipoleInRing = 1
+pulsar.DampBDipoleInRing = 0
 pulsar.Bdamping_scale = 10000
+pulsar.LimitDipoleBInit = 0   # if 1 then Bdipole is iniialization only in a smaller region
+pulsar.DipoleB_init_maxradius = 12032 # read if LimitBdipole is 1
+pulsar.use_theoreticalEB = 0
+pulsar.theory_max_rstar = 100000
+pulsar.EnforceTheoreticalEBInGrid = 1
+pulsar.AddExternalMonopoleOnly = 1
+

Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.H

@@ -16,7 +16,7 @@ public:
                    const int lev,
                    const amrex::IntVect crse_ratio,
                    int ncomp = 1);
-    virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override;    
+    virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override;
     void ComputeEdotB(amrex::MultiFab& mf_dst, int dcomp) const;
 private:
     amrex::MultiFab const * const m_Ex_src = nullptr;

Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.cpp

@@ -21,10 +21,6 @@ void
 EdotBFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const
 {
     using namespace amrex;
-    auto & warpx = WarpX::GetInstance();
-    const auto dx = warpx.Geom(m_lev).CellSizeArray();
-    const auto problo = warpx.Geom(m_lev).ProbLoArray();
-    const auto probhi = warpx.Geom(m_lev).ProbHiArray();
     const amrex::IntVect stag_dst = mf_dst.ixType().toIntVect();
 
     // convert boxarray of source MultiFab to staggering of dst Multifab
@@ -38,25 +34,21 @@ EdotBFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buff
         ComputeEdotB(mf_dst, dcomp);
     } else {
         const int ncomp = 1;
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(), 
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(),
             " all sources must have the same Distribution map");
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src-> DistributionMap() == m_Bz_src->DistributionMap(), 
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src-> DistributionMap() == m_Bz_src->DistributionMap(),
             " all sources must have the same Distribution map");
         amrex::MultiFab mf_tmp( ba_tmp, m_Ex_src->DistributionMap(), ncomp, 0);
         const int dcomp_tmp = 0;
         ComputeEdotB(mf_tmp, dcomp_tmp);
-        mf_dst.copy( mf_tmp, 0, dcomp, ncomp);
+        mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp);
     }
 }
 
 void
 EdotBFunctor::ComputeEdotB(amrex::MultiFab& mf_dst, int dcomp) const
 {
     using namespace amrex;
-    auto & warpx = WarpX::GetInstance();
-    const auto dx = warpx.Geom(m_lev).CellSizeArray();
-    const auto problo = warpx.Geom(m_lev).ProbLoArray();
-    const auto probhi = warpx.Geom(m_lev).ProbHiArray();
     const amrex::IntVect stag_Exsrc = m_Ex_src->ixType().toIntVect();
     const amrex::IntVect stag_Eysrc = m_Ey_src->ixType().toIntVect();
     const amrex::IntVect stag_Ezsrc = m_Ez_src->ixType().toIntVect();

Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.H

@@ -17,7 +17,7 @@ public:
                             const amrex::IntVect crse_ratio,
                             int vectorcomp,
                             int ncomp = 1);
-    virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override;    
+    virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override;
 
     void ComputePoyntingVector(amrex::MultiFab& mf_dst, int dcomp) const;
 private:

Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.cpp

@@ -25,10 +25,6 @@ PoyntingVectorFunctor::operator ()(amrex::MultiFab& mf_dst,
                                    int dcomp, const int /*i_buffer=0*/) const
 {
     using namespace amrex;
-    auto & warpx = WarpX::GetInstance();
-    const auto dx = warpx.Geom(m_lev).CellSizeArray();
-    const auto problo = warpx.Geom(m_lev).ProbLoArray();
-    const auto probhi = warpx.Geom(m_lev).ProbHiArray();
     const amrex::IntVect stag_dst = mf_dst.ixType().toIntVect();
 
     // convert boxarray of source MultiFab to staggering of dst Multifab
@@ -42,14 +38,14 @@ PoyntingVectorFunctor::operator ()(amrex::MultiFab& mf_dst,
         ComputePoyntingVector(mf_dst, dcomp);
     } else {
         const int ncomp = 1;
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(), 
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(),
             " all sources must have the same Distribution map");
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src->DistributionMap() == m_Bz_src->DistributionMap(), 
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src->DistributionMap() == m_Bz_src->DistributionMap(),
             " all sources must have the same Distribution map");
         amrex::MultiFab mf_tmp( ba_tmp, m_Ex_src->DistributionMap(), ncomp, 0 );
         const int dcomp_tmp = 0;
         ComputePoyntingVector(mf_tmp, dcomp_tmp);
-        mf_dst.copy( mf_tmp, 0, dcomp, ncomp);
+        mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp);
     }
 
 }
@@ -58,10 +54,6 @@ void
 PoyntingVectorFunctor::ComputePoyntingVector(amrex::MultiFab& mf_dst, int dcomp) const
 {
     using namespace amrex;
-    auto & warpx = WarpX::GetInstance();
-    const auto dx = warpx.Geom(m_lev).CellSizeArray();
-    const auto problo = warpx.Geom(m_lev).ProbLoArray();
-    const auto probhi = warpx.Geom(m_lev).ProbHiArray();
     const amrex::IntVect stag_Exsrc = m_Ex_src->ixType().toIntVect();
     const amrex::IntVect stag_Eysrc = m_Ey_src->ixType().toIntVect();
     const amrex::IntVect stag_Ezsrc = m_Ez_src->ixType().toIntVect();
@@ -90,7 +82,6 @@ PoyntingVectorFunctor::ComputePoyntingVector(amrex::MultiFab& mf_dst, int dcomp)
         cr[i] = m_crse_ratio[i];
     }
     const int vectorcomp = m_vectorcomp;
-    amrex::Real mu0_inv = 1.0_rt/PhysConst::mu0;
 
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())

Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.H

@@ -8,7 +8,7 @@ SphericalComponentFunctor final : public ComputeDiagFunctor
 {
 public:
     SphericalComponentFunctor ( const amrex::MultiFab * const mfx_src,
-                                const amrex::MultiFab * const mfy_src, 
+                                const amrex::MultiFab * const mfy_src,
                                 const amrex::MultiFab * const mfz_src,
                                 const int lev,
                                 const amrex::IntVect crse_ratio,
@@ -26,7 +26,7 @@ private:
     amrex::MultiFab const * const m_mfz_src = nullptr;
     int m_lev;
     int m_sphericalcomp;
-    int m_Efield;    
+    int m_Efield;
 };
 
 #endif

Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.cpp

@@ -12,7 +12,7 @@
 #include <AMReX.H>
 
 SphericalComponentFunctor::SphericalComponentFunctor (amrex::MultiFab const * mfx_src,
-                                                      amrex::MultiFab const * mfy_src, 
+                                                      amrex::MultiFab const * mfy_src,
                                                       amrex::MultiFab const * mfz_src,
                                                       int lev,
                                                       amrex::IntVect crse_ratio,
@@ -25,14 +25,10 @@ SphericalComponentFunctor::SphericalComponentFunctor (amrex::MultiFab const * mf
 {}
 
 
-void 
+void
 SphericalComponentFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const
 {
     using namespace amrex;
-    auto & warpx = WarpX::GetInstance();
-    const auto dx = warpx.Geom(m_lev).CellSizeArray();
-    const auto problo = warpx.Geom(m_lev).ProbLoArray();
-    const auto probhi = warpx.Geom(m_lev).ProbHiArray();
     const amrex::IntVect stag_dst = mf_dst.ixType().toIntVect();
 
     // convert boxarray of source MultiFab to staggering of dst Multifab
@@ -46,23 +42,23 @@ SphericalComponentFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const
         ComputeSphericalFieldComponent(mf_dst, dcomp);
     } else {
         const int ncomp = 1;
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_mfx_src->DistributionMap() == m_mfy_src->DistributionMap() and m_mfy_src->DistributionMap() == m_mfz_src->DistributionMap(), 
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_mfx_src->DistributionMap() == m_mfy_src->DistributionMap() and m_mfy_src->DistributionMap() == m_mfz_src->DistributionMap(),
             " all sources must have the same Distribution map");
         amrex::MultiFab mf_tmp( ba_tmp, m_mfx_src->DistributionMap(), ncomp, 0);
         const int dcomp_tmp = 0;
         ComputeSphericalFieldComponent(mf_tmp, dcomp_tmp);
-        mf_dst.copy( mf_tmp, 0, dcomp, ncomp);
+        mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp);
     }
 }
 
 void
 SphericalComponentFunctor::ComputeSphericalFieldComponent( amrex::MultiFab& mf_dst, int dcomp) const
 {
     using namespace amrex;
+#ifdef PULSAR
     auto & warpx = WarpX::GetInstance();
     const auto dx = warpx.Geom(m_lev).CellSizeArray();
     const auto problo = warpx.Geom(m_lev).ProbLoArray();
-    const auto probhi = warpx.Geom(m_lev).ProbHiArray();
     const amrex::IntVect stag_xsrc = m_mfx_src->ixType().toIntVect();
     const amrex::IntVect stag_ysrc = m_mfy_src->ixType().toIntVect();
     const amrex::IntVect stag_zsrc = m_mfz_src->ixType().toIntVect();
@@ -73,18 +69,16 @@ SphericalComponentFunctor::ComputeSphericalFieldComponent( amrex::MultiFab& mf_d
     amrex::GpuArray<int,3> sfz; // staggering of source zfield
     amrex::GpuArray<int,3> s_dst;
     amrex::GpuArray<int,3> cr;
-
+    amrex::GpuArray<amrex::Real, 3> center_star_arr;
     for (int i=0; i<AMREX_SPACEDIM; ++i) {
         sfx[i] = stag_xsrc[i];
         sfy[i] = stag_ysrc[i];
         sfz[i] = stag_zsrc[i];
         s_dst[i]  = stag_dst[i];
         cr[i] = m_crse_ratio[i];
+        center_star_arr[i] = Pulsar::m_center_star[i];
     }
     const int sphericalcomp = m_sphericalcomp;
-    amrex::Real cur_time = warpx.gett_new(0);
-    int Efield = m_Efield;
-#ifdef PULSAR
 
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
@@ -109,24 +103,25 @@ SphericalComponentFunctor::ComputeSphericalFieldComponent( amrex::MultiFab& mf_d
                 // convert to spherical coordinates
                 // compute cell coordinates
                 amrex::Real x, y, z;
-                PulsarParm::ComputeCellCoordinates(i,j,k, s_dst, problo, dx, x, y, z);
+                Pulsar::ComputeCellCoordinates(i,j,k, s_dst, problo, dx, x, y, z);
                 // convert cartesian to spherical coordinates
                 amrex::Real r, theta, phi;
-                PulsarParm::ConvertCartesianToSphericalCoord(x, y, z, problo, probhi,
-                                                             r, theta, phi);
+                Pulsar::ConvertCartesianToSphericalCoord(x, y, z, center_star_arr,
+                                                         r, theta, phi);
 
-                if (sphericalcomp == 0) { // rcomponent of field 
-                    PulsarParm::ConvertCartesianToSphericalRComponent(
+                if (sphericalcomp == 0) { // rcomponent of field
+                    Pulsar::ConvertCartesianToSphericalRComponent(
                         cc_xfield, cc_yfield, cc_zfield, theta, phi, arr_dst(i,j,k,n+dcomp));
                 } else if (sphericalcomp == 1) { // theta component of field
-                    PulsarParm::ConvertCartesianToSphericalThetaComponent(
+                    Pulsar::ConvertCartesianToSphericalThetaComponent(
                         cc_xfield, cc_yfield, cc_zfield, theta, phi, arr_dst(i,j,k,n+dcomp));
                 } else if (sphericalcomp == 2) { // phi component of field
-                    PulsarParm::ConvertCartesianToSphericalPhiComponent(
+                    Pulsar::ConvertCartesianToSphericalPhiComponent(
                         cc_xfield, cc_yfield, cc_zfield, theta, phi, arr_dst(i,j,k,n+dcomp));
                 }
             });
     }
-
+#else
+    amrex::ignore_unused(mf_dst,dcomp);
 #endif
 }

Source/Diagnostics/ReducedDiags/FieldReduction.H

@@ -67,6 +67,7 @@ private:
 
     // Type of reduction (e.g. Maximum, Minimum or Sum)
     int m_reduction_type;
+    int m_integral_type;
 
 public:
 
@@ -128,6 +129,7 @@ public:
 
         // get parser
         auto reduction_function_parser = m_parser->compile<m_nvars>();
+        int integral_type = m_integral_type;
 
         // MFIter loop to interpolate fields to cell center and perform reduction
 #ifdef AMREX_USE_OMP
@@ -190,11 +192,20 @@ public:
             amrex::ParallelDescriptor::ReduceRealSum(reduce_value);
         // If reduction operation is a sum, multiply the value by the cell volume so that the
         // result is the integral of the function over the simulation domain.
+            if (integral_type == 0) {
 #if (AMREX_SPACEDIM==2)
             reduce_value *= dx[0]*dx[1];
 #else
             reduce_value *= dx[0]*dx[1]*dx[2];
 #endif
+            } else {
+            // works only when dx = dy = dz
+#if (AMREX_SPACEDIM==2)
+            reduce_value *= dx[0];
+#else
+            reduce_value *= dx[0]*dx[1];
+#endif
+            }
         }
 
         // Fill output array

Source/Diagnostics/ReducedDiags/FieldReduction.cpp

@@ -61,6 +61,9 @@ FieldReduction::FieldReduction (std::string rd_name)
     std::string reduction_type_string;
     pp_rd_name.get("reduction_type", reduction_type_string);
     m_reduction_type = GetAlgorithmInteger (pp_rd_name, "reduction_type");
+    if (m_reduction_type == 2) {
+        m_integral_type = GetAlgorithmInteger (pp_rd_name, "integration_type");
+    }
 
     if (amrex::ParallelDescriptor::IOProcessor())
     {