Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Pulsar add obliquity #97

Open
wants to merge 13 commits into
base: NewPulsarDevelopment
Choose a base branch
from
Open

Pulsar add obliquity #97

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
945312f
add ndens and magnetization mf, visualize, add particles if local si…
RevathiJambunathan Apr 7, 2022
74a9bba
add plasma if local magnetization is higher that threshold
RevathiJambunathan Apr 7, 2022
2a0c3e2
fix restart by adding sum_injection_rate to IO
RevathiJambunathan Apr 28, 2022
440f636
typo
RevathiJambunathan Apr 28, 2022
ab108de
for magnetization use cell-centered B, and compute sum of injection flag
RevathiJambunathan May 9, 2022
28b6849
print injected cell details
RevathiJambunathan May 9, 2022
6ac7fcf
print only data for injected cells
RevathiJambunathan May 10, 2022
c30e2cc
add Barrier and new var for sigma0
RevathiJambunathan May 11, 2022
dbac046
remove vismf write
RevathiJambunathan May 12, 2022
7a7170b
add comments, cleanup, and add user-defined var for ubound of rel_dif…
RevathiJambunathan May 16, 2022
2bd3a5b
var for obliquity
RevathiJambunathan May 16, 2022
37e0d76
add oblique angle Chi based on Eq 2.7, 2.8, 2.9 in Jeromi Petris 2016…
RevathiJambunathan May 17, 2022
42bac1f
fix tab
RevathiJambunathan May 19, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions Source/Evolve/WarpXEvolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down
10 changes: 6 additions & 4 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Open in desktop
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
140 changes: 101 additions & 39 deletions Source/Particles/PulsarParameters.H
View file
Open in desktop
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 Expand Up @@ -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:
};

Expand Down
Loading