updated pgen

Author: StaticObserver

pgens/1d_polar_cap/pgen.hpp

@@ -19,23 +19,29 @@ namespace user {
 
   template <Dimension D>
   struct InitFields {
-    InitFields(real_t b0_, real_t rho_GJ_, real_t l_atm_, real_t ds_) 
+    InitFields(real_t b0_, real_t coeff_, real_t xsurf_, real_t ds_) 
           : b0 { b0_ }
-          , rho_GJ { rho_GJ_ } 
-          , l_atm { l_atm_ }
+          , coeff { coeff_ } 
+          , xsurf { xsurf_ }
           , ds { ds_ }{}
 
     Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
       return b0;
     }
 
     Inline auto ex1(const coord_t<D>& x_Ph) const -> real_t {
-        return -rho_GJ * (x_Ph[0] + 0.03 * ds * math::log(0.01 + math::exp(l_atm + 0.8 * ds - x_Ph[0]) / 0.03 / ds));
-        // return -rho_GJ * (x_Ph[0] - l_atm);
+          if (x_Ph[0] <  xsurf ){
+              return ZERO;
+          }else{
+               //return -coeff * (x_Ph[0] - xsurf
+               //                      + 0.03 * ds * math::log(0.01 + math::exp((xsurf + 1.0 * ds - x_Ph[0]) / 0.03 / ds))
+                 //                    - 0.03 * ds * math::log(0.01 + math::exp(1.0 / 0.03)));
+              return -coeff * (x_Ph[0] - xsurf - ds);
+          }
     }
 
   private:
-    const real_t b0, rho_GJ, l_atm, ds;
+    const real_t b0, coeff, xsurf, ds;
   };
 
   template <Dimension D>
@@ -57,8 +63,8 @@ namespace user {
 
   template <Dimension D>
   struct DriveFields : public InitFields<D> {
-    DriveFields(real_t time, real_t b0_, real_t l_atm_, real_t ds_)
-      : InitFields<D> { b0_, ZERO, l_atm_, ds_ } {}
+    DriveFields(real_t time, real_t b0_, real_t xsurf_, real_t ds_)
+      : InitFields<D> { b0_, ZERO, xsurf_, ds_ } {}
 
     using InitFields<D>::bx1;
 
@@ -71,17 +77,22 @@ namespace user {
 
   template <SimEngine::type S, class M>
     struct TargetDensityProfile : public arch::SpatialDistribution<S, M> {
-        const real_t nmax, height, xsurf;
+        const real_t nmax, height, xsurf, ds;
 
-        TargetDensityProfile(const M& metric, real_t nmax, real_t height, real_t xsurf)
+        TargetDensityProfile(const M& metric, real_t nmax, real_t height, real_t xsurf, real_t ds)
           : arch::SpatialDistribution<S, M>(metric)
           , nmax { nmax }
           , height { height }
-          , xsurf { xsurf } {}
+          , xsurf { xsurf }
+          , ds { ds } {}
 
         Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t {
-              return nmax * math::exp(-(x_Ph[0] - xsurf) / height);
-          }
+             if (x_Ph[0] < xsurf){
+                  return ZERO;
+             }else{
+                  return nmax * math::exp(-(x_Ph[0] - xsurf) / height);
+             }
+           }
       }; // TargetDensityProfile
 
     template <SimEngine::type S, class M>
@@ -94,7 +105,11 @@ namespace user {
           , ds { ds_ } {}  
 
         Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t {
-              return ONE - 0.01 / (0.01 + math::exp(-(x_Ph[0] - xsurf - 0.8 * ds) / 0.03 / ds));
+               if ( x_Ph[0] < xsurf ){
+                  return ZERO;
+               }else{
+                   return ONE - 0.01 / (0.01 + math::exp(-(x_Ph[0] - xsurf - 1.0 * ds) / 0.03 / ds));
+               }
           }
       }; // ExtraCharge
 
@@ -137,7 +152,7 @@ namespace user {
     const real_t  b0, Omega, skindepth0, larmor0;
     const real_t  temp;
     const real_t  j0;
-    const real_t  l_atm, ds;
+    const real_t  xsurf, ds;
     InitFields<D> init_flds;
     MagnetosphericCurrent<D> ext_current;
 
@@ -150,23 +165,24 @@ namespace user {
       , skindepth0 { p.template get<real_t>("scales.skindepth0") }
       , larmor0 { p.template get<real_t>("scales.larmor0") }
       , temp { p.template get<real_t>("setup.temp") }
-      , j0 { p.template get<real_t>("setup.j0") }
+      , j0 { p.template get<real_t>("setup.j0") / p.template get<real_t>("scales.V0")}
       , ds { p.template get<real_t>("grid.boundaries.atmosphere.ds") }
-      , l_atm([&m, this]() {
+      , xsurf([&m, this]() {
           const auto min_buff = params.template get<unsigned short>("algorithms.current_filters") + 2;
           const auto buffer_ncells = min_buff > 5 ? min_buff : 5;
           return m.mesh().metric.template convert<1, Crd::Cd, Crd::Ph>(static_cast<real_t>(buffer_ncells));
         }())
+      //, l_atm { ZERO } 
       // , init_flds(b0, TWO * FOUR * constant::PI * b0 * Omega * SQR(skindepth0) / larmor0, l_atm, ds)
-      , init_flds(b0, ONE, l_atm, ds)
+      , init_flds(b0, larmor0 / SQR(skindepth0), xsurf, ds)
       , ext_current(j0)
     {}
 
     inline PGen() {}
 
 
     auto AtmFields(real_t time) const -> DriveFields<D> {
-      return DriveFields<D> { time, b0, l_atm, ds};
+      return DriveFields<D> { time, b0, xsurf, ds};
     }
 
     auto MatchFields(real_t) const -> MFields<D> {
@@ -193,7 +209,8 @@ namespace user {
           local_domain.mesh.metric,
           params.template get<real_t>("grid.boundaries.atmosphere.density"),
           params.template get<real_t>("grid.boundaries.atmosphere.height"),
-          l_atm);
+          xsurf,
+          ds);
       const auto injector = arch::NonUniformInjector<S, M, arch::Maxwellian, TargetDensityProfile>(
         energy_dist,
         spatial_dist,
@@ -207,19 +224,19 @@ namespace user {
 
       const auto extra_charge = ExtraCharge<S, M>(
         local_domain.mesh.metric,
-        l_atm,
+        xsurf,
         ds
       );
       const auto injector_extra_charge = arch::NonUniformInjector<S, M, arch::Maxwellian, ExtraCharge>(
         energy_dist,
         extra_charge,
         { 2, 2 }
       );
-      arch::InjectNonUniform<S, M, decltype(injector_extra_charge)>(
-        params,
-        local_domain,
-        injector_extra_charge,
-        TWO);
+      //arch::InjectNonUniform<S, M, decltype(injector_extra_charge)>(
+       // params,
+       // local_domain,
+       // injector_extra_charge,
+       // TWO);
     }
   }; // PGen