diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index c2d00b5de66..95eb288a3f4 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -354,6 +354,12 @@ WarpX::Evolve (int numsteps) } } +#ifdef PULSAR + m_pulsar->ComputePlasmaNumberDensity(); + m_pulsar->ComputePlasmaMagnetization(); +#endif + + if (sort_intervals.contains(step+1)) { if (verbose) { amrex::Print() << Utils::TextMsg::Info("re-sorting particles"); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index cfba3cb7541..e9c2e5c93f6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2320,6 +2320,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, amrex::Real Bstar_data = Pulsar::m_B_star; amrex::Real Rstar_data = Pulsar::m_R_star; amrex::Real dRstar_data = Pulsar::m_dR_star; + amrex::Real chi_data = Pulsar::m_Chi; amrex::Real corotatingE_maxradius_data = Pulsar::m_corotatingE_maxradius; int E_external_monopole_data = Pulsar::m_do_E_external_monopole; int AddExternalMonopoleOnly = Pulsar::m_AddExternalMonopoleOnly; @@ -2425,7 +2426,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, amrex::ParticleReal r_p, theta_p, phi_p; Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr, r_p, theta_p, phi_p); - Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p, + Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p, chi_data, AddExternalMonopoleOnly, AddPulsarVacuumEFields, AddBDipoleOutsideRstar, @@ -2436,7 +2437,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, Bstar_data, Rstar_data, dRstar_data, Exp, Eyp, Ezp, Bxp, Byp, Bzp); if (use_theoreticalEB == 1) { - Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p, + Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p, chi_data, theory_max_rstar, corotatingE_maxradius_data, E_external_monopole_data, @@ -2747,6 +2748,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, amrex::Real Bstar_data = Pulsar::m_B_star; amrex::Real Rstar_data = Pulsar::m_R_star; amrex::Real dRstar_data = Pulsar::m_dR_star; + amrex::Real chi_data = Pulsar::m_Chi; amrex::Real corotatingE_maxradius_data = Pulsar::m_corotatingE_maxradius; int E_external_monopole_data = Pulsar::m_do_E_external_monopole; int AddExternalMonopoleOnly = Pulsar::m_AddExternalMonopoleOnly; @@ -2878,7 +2880,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, amrex::ParticleReal r_p, theta_p, phi_p; Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr, r_p, theta_p, phi_p); - Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p, + Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p, chi_data, AddExternalMonopoleOnly, AddPulsarVacuumEFields, AddBDipoleOutsideRstar, @@ -2889,7 +2891,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, Bstar_data, Rstar_data, dRstar_data, Exp, Eyp, Ezp, Bxp, Byp, Bzp); if (use_theoreticalEB == 1) { - Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p, + Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p, chi_data, theory_max_rstar, corotatingE_maxradius_data, E_external_monopole_data, diff --git a/Source/Particles/PulsarParameters.H b/Source/Particles/PulsarParameters.H index f3b24420972..a81d63c242d 100644 --- a/Source/Particles/PulsarParameters.H +++ b/Source/Particles/PulsarParameters.H @@ -58,7 +58,8 @@ public: AMREX_GPU_HOST_DEVICE AMREX_INLINE static void CorotatingEfieldSpherical (amrex::Real const r, amrex::Real const theta, - amrex::Real const phi, amrex::Real const time, + amrex::Real const phi, amrex::Real const chi, + amrex::Real const time, amrex::Real const omega_star_data, amrex::Real const ramp_omega_time_data, amrex::Real const Bstar, @@ -66,10 +67,18 @@ public: amrex::Real const dRstar, amrex::Real &Er, amrex::Real &Etheta, amrex::Real &Ephi) { - amrex::ignore_unused(phi); amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data); + // Polar angle amrex::Real c_theta = std::cos(theta); amrex::Real s_theta = std::sin(theta); + // Oblique angle + amrex::Real c_chi = std::cos(chi); + amrex::Real s_chi = std::sin(chi); + // Instantaneous phase, psi = phi - Omega*t (phi is the azimuthal angle) + // Refer to Pg. 6 of Jeromi Petri's 2016 paper + amrex::Real psi = phi - omega*time; + amrex::Real c_psi = std::cos(psi); + // Ratio : Rstar/r amrex::Real r_ratio; if (r > 0) { r_ratio = Rstar/r; @@ -77,15 +86,21 @@ public: r_ratio = Rstar/(dRstar*0.5); } amrex::Real r2 = r_ratio * r_ratio; - // Michel and Li -- eq 14 , 15 - Er = Bstar * omega * r2 * Rstar * s_theta * s_theta; - Etheta = -Bstar * omega * r2 * Rstar * 2.0 * s_theta * c_theta; - Ephi = 0.0; // aligned magnetic and rotation axis + // Corotating electric field inside the pulsar that corresponds to dipole magnetic field + // Eq 2.9 (a-c) Jeromi Petri 2016 + Er = Bstar * omega * r2 * Rstar * + ( c_chi * s_theta * s_theta + - s_chi * c_theta * s_theta * c_psi); + Etheta = -Bstar * omega * r2 * Rstar * + ( c_chi * 2._rt * s_theta * c_theta + + s_chi * 2._rt * s_theta * s_theta * c_psi); + Ephi = 0._rt; } AMREX_GPU_HOST_DEVICE AMREX_INLINE static void ExternalEFieldSpherical (amrex::Real const r, amrex::Real const theta, - amrex::Real const phi, amrex::Real const time, + amrex::Real const phi, amrex::Real const chi, + amrex::Real const time, amrex::Real const omega_star_data, amrex::Real const ramp_omega_time_data, amrex::Real const Bstar, amrex::Real const Rstar, @@ -94,38 +109,59 @@ public: int const ApplyCorotatingEField, amrex::Real &Er, amrex::Real &Etheta, amrex::Real &Ephi) { - amrex::ignore_unused(phi, corotatingE_maxradius); + amrex::ignore_unused(corotatingE_maxradius); amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data); + // Polar angle amrex::Real c_theta = std::cos(theta); amrex::Real s_theta = std::sin(theta); + // Oblique angle + amrex::Real c_chi = std::cos(chi); + amrex::Real s_chi = std::sin(chi); + // Instantaneous phase, psi = phi - Omega*t (phi is the azimuthal angle) + // Refer to Pg. 6 of Jeromi Petri's 2016 paper + amrex::Real psi = phi - omega*time; + amrex::Real c_psi = std::cos(psi); + amrex::Real s_psi = std::sin(psi); + amrex::Real r_ratio = Rstar/r; // inside pulsar - //if (r <= corotatingE_maxradius) { if (ApplyCorotatingEField == 1) { amrex::Real r2 = r_ratio * r_ratio; - // Michel and Li -- eq 14 , 15 - Er = Bstar * omega * r2 * Rstar * s_theta * s_theta; - Etheta = -Bstar * omega * r2 * Rstar * 2.0 * s_theta * c_theta; - Ephi = 0.0; // aligned magnetic and rotation axis + // Corotating electric field inside the pulsar that corresponds to dipole magnetic field + // Eq 2.9 (a-c) Jeromi Petri 2016 + Er = Bstar * omega * r2 * Rstar * + ( c_chi * s_theta * s_theta + - s_chi * c_theta * s_theta * c_psi); + Etheta = -Bstar * omega * r2 * Rstar * + ( c_chi * 2._rt * s_theta * c_theta + + s_chi * 2._rt * s_theta * s_theta * c_psi); + Ephi = 0._rt; } // outside pulsar - //if (r > corotatingE_maxradius) { if (ApplyCorotatingEField == 0) { - amrex::Real r4 = r_ratio*r_ratio*r_ratio*r_ratio; - // Taking derivative of phi given in eq 30 of Michel and Li - Er = Bstar * omega * Rstar * r4 * (1.0-3.0*c_theta*c_theta); + amrex::Real r2 = r_ratio * r_ratio; + amrex::Real r4 = r2 * r2; + // Equations 2.8(a) - 2.8(c) + Er = Bstar * omega * Rstar * r4 * + ( c_chi * (1._rt - 3._rt * c_theta * c_theta) + - 3._rt * s_chi * c_theta * s_theta * c_psi ); if (Eexternal_monopole == 1) { - Er += (2.0/3.0) * omega * Bstar * Rstar * r_ratio * r_ratio; + // Monopole term from central charge = Qc/(4*pi*eps0*r^2) + // For dipolar magnetic field, Qc = 8pi/3 * eps0*Omega*Bstar*Rstar^3*c_chi + Er += (2._rt/3._rt) * omega * Bstar * Rstar * r2 * c_chi; } - Etheta = (-1.0) * Bstar * omega * Rstar * r4 * (2.0*s_theta*c_theta); - Ephi = 0.0; + Etheta = Bstar * omega * Rstar * + ( r2 * s_chi * (r2*std::cos(2._rt*theta) - 1._rt) * c_psi + - r4 * c_chi * std::sin(2._rt*theta) ); + Ephi = Bstar * omega * Rstar * r2 * (1._rt - r2) * s_chi * c_theta * s_psi; } } AMREX_GPU_HOST_DEVICE AMREX_INLINE static void ExternalEMonopoleSpherical (amrex::Real const r, amrex::Real const theta, - amrex::Real const phi, amrex::Real const time, + amrex::Real const phi, amrex::Real const chi, + amrex::Real const time, amrex::Real const omega_star_data, amrex::Real const ramp_omega_time_data, amrex::Real const Bstar, amrex::Real const Rstar, @@ -134,8 +170,10 @@ public: amrex::ignore_unused(phi, theta); amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data); amrex::Real r_ratio = Rstar/r; + amrex::Real r2 = r_ratio * r_ratio; + amrex::Real c_chi = std::cos(chi); - Er = (2.0/3.0) * omega * Bstar * Rstar * r_ratio * r_ratio; + Er = (2._rt/3._rt) * omega * Bstar * Rstar * r2 * c_chi; Etheta = 0.; Ephi = 0.; @@ -143,14 +181,27 @@ public: AMREX_GPU_HOST_DEVICE AMREX_INLINE static void ExternalBFieldSpherical (amrex::Real const r, amrex::Real const theta, - amrex::Real const phi, amrex::Real const time, + amrex::Real const phi, amrex::Real const chi, + amrex::Real const time, + amrex::Real const omega_star_data, + amrex::Real const ramp_omega_time_data, amrex::Real const Bstar, amrex::Real const Rstar, amrex::Real const dRstar, amrex::Real &Br, amrex::Real &Btheta, amrex::Real &Bphi) { - amrex::ignore_unused(phi, time); + // Polar angle amrex::Real c_theta = std::cos(theta); amrex::Real s_theta = std::sin(theta); + // Oblique angle + amrex::Real c_chi = std::cos(chi); + amrex::Real s_chi = std::sin(chi); + // Instantaneous phase, psi = phi - Omega*t (phi is the azimuthal angle) + // Refer to Pg. 6 of Jeromi Petri's 2016 paper + amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data); + amrex::Real psi = phi - omega*time; + amrex::Real c_psi = std::cos(psi); + amrex::Real s_psi = std::sin(psi); + amrex::Real r_ratio; if (r > 0) { r_ratio = Rstar/r; @@ -158,15 +209,16 @@ public: r_ratio = Rstar/(dRstar*0.5); } amrex::Real r3 = r_ratio*r_ratio*r_ratio; - // Michel and Li -- eq 14 and 15 from michel and li - // dipole B field inside and outside the pulsar - Br = 2.0*Bstar*r3*c_theta; - Btheta = Bstar*r3*s_theta; - Bphi = 0.0; + // The full dipole magnetic field as f(theta, phi, omega, t, and Chi) + // Eqs 2.7(a) - 2.7(c) from Jeromi Petri 2016 paper + Br = 2._rt * Bstar * r3 * ( c_chi * c_theta + s_chi * s_theta * c_psi); + Btheta = Bstar * r3 * ( c_chi * s_theta - s_chi * c_theta * c_psi); + Bphi = Bstar * r3 * s_chi * s_psi; } AMREX_GPU_HOST_DEVICE AMREX_INLINE - static void getExternalEBOnParticle( amrex::Real const r, amrex::Real const theta, amrex::Real const phi, + static void getExternalEBOnParticle( amrex::Real const r, amrex::Real const theta, + amrex::Real const phi, amrex::Real const chi, const int AddExternalMonopoleOnly, const int AddPulsarVacuumEFields, const int AddBDipoleOutsideRstar, const int AddPulsarVacuumBFields, amrex::Real const corotatingE_maxradius, @@ -184,7 +236,7 @@ public: amrex::Real Er = 0.0_rt; amrex::Real Etheta = 0.0_rt; amrex::Real Ephi = 0.0_rt; - ExternalEMonopoleSpherical(r, theta, phi, cur_time, omega_star, + ExternalEMonopoleSpherical(r, theta, phi, chi, cur_time, omega_star, ramp_omega_time, Bstar, Rstar, Er, Etheta, Ephi); amrex::Real Ex_monopole = 0._rt; @@ -206,7 +258,7 @@ public: amrex::Real Ephi = 0.0_rt; int ApplyCorotatingEField = 0; if (r <= corotatingE_maxradius) ApplyCorotatingEField = 1; - ExternalEFieldSpherical(r, theta, phi, cur_time, omega_star, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star, ramp_omega_time, Bstar, Rstar, corotatingE_maxradius, E_external_monopole, @@ -229,8 +281,10 @@ public: amrex::Real Br = 0._rt; amrex::Real Btheta = 0._rt; amrex::Real Bphi = 0._rt; - ExternalBFieldSpherical (r, theta, phi, cur_time, Bstar, Rstar, dR_star, - Br, Btheta, Bphi); + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star, ramp_omega_time, + Bstar, Rstar, dR_star, + Br, Btheta, Bphi); amrex::Real Bx_dipole = 0._rt; amrex::Real By_dipole = 0._rt; amrex::Real Bz_dipole = 0._rt; @@ -246,8 +300,10 @@ public: amrex::Real Br = 0._rt; amrex::Real Btheta = 0._rt; amrex::Real Bphi = 0._rt; - ExternalBFieldSpherical (r, theta, phi, cur_time, Bstar, Rstar, dR_star, - Br, Btheta, Bphi); + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star, ramp_omega_time, + Bstar, Rstar, dR_star, + Br, Btheta, Bphi); amrex::Real Bx_dipole = 0._rt; amrex::Real By_dipole = 0._rt; amrex::Real Bz_dipole = 0._rt; @@ -262,7 +318,8 @@ public: AMREX_GPU_HOST_DEVICE AMREX_INLINE - static void EnforceTheoreticalEBOnParticle( amrex::Real const r, amrex::Real const theta, amrex::Real const phi, + static void EnforceTheoreticalEBOnParticle( amrex::Real const r, amrex::Real const theta, + amrex::Real const phi, amrex::Real const chi, amrex::Real const theory_max_radius, amrex::Real const corotatingE_maxradius, const int E_external_monopole, @@ -285,7 +342,7 @@ public: amrex::Real Ex_theory = 0._rt; amrex::Real Ey_theory = 0._rt; amrex::Real Ez_theory = 0._rt; - ExternalEFieldSpherical(r, theta, phi, cur_time, omega_star, ramp_omega_time, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star, ramp_omega_time, Bstar, Rstar, corotatingE_maxradius, E_external_monopole, ApplyCorotatingEField, Er, Etheta, Ephi); // Compute theoretical Bfield in spherical coordinates @@ -295,7 +352,8 @@ public: amrex::Real Bx_theory = 0._rt; amrex::Real By_theory = 0._rt; amrex::Real Bz_theory = 0._rt; - ExternalBFieldSpherical(r, theta, phi, cur_time, Bstar, Rstar, dR_star, Br, Btheta, Bphi); + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, omega_star, ramp_omega_time, + Bstar, Rstar, dR_star, Br, Btheta, Bphi); // Convert Efield from spherical to Cartesian ConvertSphericalToCartesianXComponent(Er, Etheta, Ephi, r, theta, phi, Ex_theory); ConvertSphericalToCartesianYComponent(Er, Etheta, Ephi, r, theta, phi, Ey_theory); @@ -529,6 +587,10 @@ public: ** Upper bound specified must be less than 1, and default is 0.1 (10%) */ static amrex::Real m_ubound_reldiff_sigma0; + /** Oblique angle between magnetic field axis and pulsar rotation axis. + ** (Read in degrees and converted to radians) + */ + static amrex::Real m_Chi; private: }; diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp index 8d1d7a9aa4e..d1127033676 100644 --- a/Source/Particles/PulsarParameters.cpp +++ b/Source/Particles/PulsarParameters.cpp @@ -91,6 +91,7 @@ int Pulsar::m_print_injected_celldata; int Pulsar::m_print_celldata_starttime; amrex::Real Pulsar::m_lbound_ndens_magnetization = 1.e-16; amrex::Real Pulsar::m_ubound_reldiff_sigma0 = 0.1; +amrex::Real Pulsar::m_Chi; Pulsar::Pulsar () @@ -100,6 +101,7 @@ Pulsar::Pulsar () void Pulsar::ReadParameters () { + amrex::Print() << " pulsar read data \n"; amrex::ParmParse pp("pulsar"); pp.query("pulsarType",m_pulsar_type); @@ -250,6 +252,11 @@ Pulsar::ReadParameters () { if (m_sigma_tune_method == "relative_difference") { pp.query("upperBound_reldiff_sigma0", m_ubound_reldiff_sigma0); } + // Obliquity of the pulsar + pp.get("Chi", m_Chi); + // convert degrees to radians + m_Chi = m_Chi * MathConst::pi / 180._rt; + amrex::Print() << "Oblique angle between B-axis and Omega-axis " << m_Chi << " radians \n"; } @@ -309,6 +316,7 @@ Pulsar::InitDataAtRestart () void Pulsar::InitData () { + amrex::Print() << " pulsar init data \n"; auto & warpx = WarpX::GetInstance(); const int nlevs_max = warpx.finestLevel() + 1; //allocate number density multifab @@ -458,6 +466,7 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult z_IndexType[idim] = z_nodal_flag[idim]; center_star_arr[idim] = m_center_star[idim]; } + amrex::Real chi = m_Chi; amrex::Real omega_star_data = m_omega_star; amrex::Real ramp_omega_time_data = m_omega_ramp_time; amrex::Real Bstar_data = m_B_star; @@ -497,12 +506,16 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult if (init_Bfield == 1) { if (LimitDipoleBInit_data == 1) { if (r < DipoleB_init_maxradius_data) { - ExternalBFieldSpherical( r, theta, phi, cur_time, + ExternalBFieldSpherical( r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } } else { - ExternalBFieldSpherical( r, theta, phi, cur_time, + ExternalBFieldSpherical( r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } @@ -512,7 +525,7 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult if (conductor(i,j,k)==1 and conductor(i+1,j,k)==1 ) { // Edge-centered Efield is fully immersed in conductor // Therefore, apply corotating electric field E = v X B - CorotatingEfieldSpherical(r, theta, phi, cur_time, + CorotatingEfieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, @@ -523,7 +536,7 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult // Apply corotating electric field E = v X B, if edge-centered Efield is fully immersed // in conductor if (conductor(i,j,k)==1 and conductor(i+1,j,k)==1) ApplyCorotatingEField = 1; - ExternalEFieldSpherical(r, theta, phi, cur_time, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, @@ -550,12 +563,16 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult if (init_Bfield == 1) { if (LimitDipoleBInit_data == 1) { if (r < DipoleB_init_maxradius_data) { - ExternalBFieldSpherical( r, theta, phi, cur_time, + ExternalBFieldSpherical( r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } } else { - ExternalBFieldSpherical( r, theta, phi, cur_time, + ExternalBFieldSpherical( r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } @@ -565,7 +582,7 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult if (conductor(i,j,k)==1 and conductor(i,j+1,k)==1 ) { // Edge-centered Efield is fully immersed in conductor // Therefore, apply corotating electric field E = v X B - CorotatingEfieldSpherical(r, theta, phi, cur_time, + CorotatingEfieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, @@ -576,7 +593,7 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult // Apply corotating electric field E = v X B, if edge-centered Efield is fully immersed // in conductor if (conductor(i,j,k)==1 and conductor(i,j+1,k)==1) ApplyCorotatingEField = 1; - ExternalEFieldSpherical(r, theta, phi, cur_time, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, @@ -603,12 +620,16 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult if (init_Bfield == 1) { if (LimitDipoleBInit_data == 1) { if (r < DipoleB_init_maxradius_data) { - ExternalBFieldSpherical( r, theta, phi, cur_time, + ExternalBFieldSpherical( r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } } else { - ExternalBFieldSpherical( r, theta, phi, cur_time, + ExternalBFieldSpherical( r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } @@ -619,7 +640,7 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult if (conductor(i,j,k)==1 and conductor(i,j,k+1)==1 ) { // Edge-centered Efield is fully immersed in conductor // Therefore, apply corotating electric field E = v X B - CorotatingEfieldSpherical(r, theta, phi, cur_time, + CorotatingEfieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, @@ -630,7 +651,7 @@ Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::Mult // Apply corotating electric field E = v X B, if edge-centered Efield is fully immersed // in conductor if (conductor(i,j,k)==1 and conductor(i,j,k+1)==1) ApplyCorotatingEField = 1; - ExternalEFieldSpherical(r, theta, phi, cur_time, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, @@ -668,6 +689,7 @@ Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, z_IndexType[idim] = z_nodal_flag[idim]; center_star_arr[idim] = m_center_star[idim]; } + amrex::Real chi = m_Chi; amrex::Real omega_star_data = m_omega_star; amrex::Real ramp_omega_time_data = m_omega_ramp_time; amrex::Real Bstar_data = m_B_star; @@ -702,7 +724,7 @@ Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, r, theta, phi); if (EnforceTheoreticalEBInGrid_data == 0) { if (conductor(i,j,k)==1 and conductor(i+1,j,k)==1) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, + CorotatingEfieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, @@ -713,7 +735,7 @@ Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, } else { int ApplyCorotatingEField = 0; if (conductor(i,j,k)==1 and conductor(i+1,j,k)==1) ApplyCorotatingEField = 1; - ExternalEFieldSpherical(r, theta, phi, cur_time, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, @@ -736,7 +758,7 @@ Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, r, theta, phi); if (EnforceTheoreticalEBInGrid_data == 0) { if (conductor(i,j,k)==1 and conductor(i,j+1,k)==1) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, + CorotatingEfieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, @@ -747,7 +769,7 @@ Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, } else { int ApplyCorotatingEField = 0; if (conductor(i,j,k)==1 and conductor(i,j+1,k)==1) ApplyCorotatingEField = 1; - ExternalEFieldSpherical(r, theta, phi, cur_time, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, @@ -771,7 +793,7 @@ Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, r, theta, phi); if (EnforceTheoreticalEBInGrid_data == 0) { if (conductor(i,j,k)==1 and conductor(i,j,k+1)==1) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, + CorotatingEfieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, @@ -782,7 +804,7 @@ Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, } else { int ApplyCorotatingEField = 0; if (conductor(i,j,k)==1 and conductor(i,j,k+1)==1) ApplyCorotatingEField = 1; - ExternalEFieldSpherical(r, theta, phi, cur_time, + ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star_data, ramp_omega_time_data, Bstar_data, Rstar_data, @@ -821,6 +843,9 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> z_IndexType[idim] = z_nodal_flag[idim]; center_star_arr[idim] = m_center_star[idim]; } + amrex::Real chi = m_Chi; + amrex::Real omega_star_data = m_omega_star; + amrex::Real ramp_omega_time_data = m_omega_ramp_time; amrex::Real Bstar_data = m_B_star; amrex::Real Rstar_data = m_R_star; amrex::Real dRstar_data = m_dR_star; @@ -858,7 +883,9 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> if (conductor(i,j,k)==1 and conductor(i,j+1,k)==1 and conductor(i,j,k+1)==1 and conductor(i,j+1,k+1)==1) { // Enforce magnetic field boundary condition (dipole field) if face-centered // Bfield location is completely embedded in the conductor - ExternalBFieldSpherical(r, theta, phi, cur_time, + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, @@ -871,9 +898,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> // to radius where corotating electric field is enforced. if (DampBDipoleInRing_data == 1) { // from inner ring to outer ring - ExternalBFieldSpherical(r, theta, phi, cur_time, - Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); amrex::Real Bx_dipole; ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, r, theta, phi, Bx_dipole); @@ -891,9 +920,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> } } else { // Theoretical Bfield enforced everywhere on the grid - ExternalBFieldSpherical(r, theta, phi, cur_time, + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + Fr, Ftheta, Fphi); ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, r, theta, phi, Bx_arr(i,j,k)); } @@ -913,9 +944,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> if (conductor(i,j,k)==1 and conductor(i+1,j,k)==1 and conductor(i,j,k+1)==1 and conductor(i+1,j,k+1)==1) { // Enforce magnetic field boundary condition (dipole field) if face-centered // Bfield location is completely embedded in the conductor - ExternalBFieldSpherical(r, theta, phi, cur_time, + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + Fr, Ftheta, Fphi); ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, r, theta, phi, By_arr(i,j,k)); } @@ -925,9 +958,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> // to radius where corotating electric field is enforced. if (DampBDipoleInRing_data == 1) { // from inner ring to outer ring - ExternalBFieldSpherical(r, theta, phi, cur_time, + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + Fr, Ftheta, Fphi); amrex::Real By_dipole; ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, r, theta, phi, By_dipole); @@ -940,9 +975,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> } } else { // Theoretical Bfield enforced everywhere on the grid - ExternalBFieldSpherical(r, theta, phi, cur_time, + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + Fr, Ftheta, Fphi); ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, r, theta, phi, By_arr(i,j,k)); } @@ -962,9 +999,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> if (conductor(i,j,k)==1 and conductor(i+1,j,k)==1 and conductor(i,j+1,k)==1 and conductor(i+1,j+1,k)==1) { // Enforce magnetic field boundary condition (dipole field) if face-centered // Bfield location is completely embedded in the conductor - ExternalBFieldSpherical(r, theta, phi, cur_time, - Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, r, theta, phi, Bz_arr(i,j,k)); } @@ -974,9 +1013,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> // to radius where corotating electric field is enforced. if (DampBDipoleInRing_data == 1) { // from inner ring to outer ring - ExternalBFieldSpherical(r, theta, phi, cur_time, - Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); amrex::Real Bz_dipole; ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, r, theta, phi, Bz_dipole); @@ -989,9 +1030,11 @@ Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> } } else { // Theoretical Bfield enforced everywhere on the grid - ExternalBFieldSpherical(r, theta, phi, cur_time, - Bstar_data, Rstar_data, dRstar_data, - Fr, Ftheta, Fphi); + ExternalBFieldSpherical(r, theta, phi, chi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, r, theta, phi, Bz_arr(i,j,k)); }