Skip to content

Commit

Permalink
add oblique angle Chi based on Eq 2.7, 2.8, 2.9 in Jeromi Petris 2016…
Browse files Browse the repository at this point in the history 
… article
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed May 20, 2022
1 parent 2bd3a5b commit 37e0d76
Show file tree
Hide file tree
Showing 3 changed files with 177 additions and 82 deletions.
10 changes: 6 additions & 4 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down
136 changes: 97 additions & 39 deletions Source/Particles/PulsarParameters.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,34 +58,49 @@ 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,
amrex::Real const Rstar,
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;
} else {
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,
Expand All @@ -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,
Expand All @@ -134,39 +170,55 @@ 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.;

}

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;
} else {
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,
Expand All @@ -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;
Expand All @@ -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,
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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,
Expand All @@ -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
Expand All @@ -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);
Expand Down
Loading

0 comments on commit 37e0d76

Please sign in to comment.