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

Jump to bottom

Minerbo initial conditions #59

Merged
merged 6 commits into from
Apr 27, 2021
Merged

Minerbo initial conditions #59

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
18d64ae
added vscode and generated_files directories to the gitignore file
Mar 23, 2021
8664070
added Minerbo closure initial conditions
Mar 23, 2021
87c00ad
added references for the minerbo equations
Mar 23, 2021
c5ad27b
added a tad more explanatory text
Mar 23, 2021
2b025a5
improved st5 code comments
Mar 23, 2021
8d428dd
fixed coefficient in front of minerbo closure to make it have an expe…
Apr 2, 2021
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
4 changes: 3 additions & 1 deletion .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,4 +37,6 @@ Scripts/symbolic_hermitians/__pycache__
.cproject
.project
.pydevproject
.settings
.settings
.vscode
Source/generated_files
150 changes: 150 additions & 0 deletions Source/FlavoredNeutrinoContainerInit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,49 @@ Gpu::ManagedVector<GpuArray<Real,3> > uniform_sphere_xyz(int nphi_at_equator){
return xyz;
}

// residual for the root finder
// Z needs to be bigger if residual is positive
dwillcox marked this conversation as resolved.
Show resolved Hide resolved
// Minerbo (1978) (unfortunately under Elsevier paywall)
// Can also see Richers (2020) https://ui.adsabs.harvard.edu/abs/2020PhRvD.102h3017R
// Eq.41 (where a is Z), but in the non-degenerate limit
// k->0, eta->0, N->Z/(4pi sinh(Z)) (just to make it integrate to 1)
// minerbo_residual is the "f" equation between eq.42 and 43
Real minerbo_residual(const Real fluxfac, const Real Z){
return fluxfac - 1.0/std::tanh(Z) + 1.0 / Z;
}
Real minerbo_residual_derivative(const Real fluxfac, const Real Z){
return 1.0/(std::sinh(Z)*std::sinh(Z)) - 1.0/(Z*Z);
}
Real minerbo_Z(const Real fluxfac){
// hard-code in these parameters because they are not
// really very important...
Real maxresidual = 1e-6;
Real maxcount = 20;
Real minfluxfac = 1e-3;

// set the initial conditions
Real Z = 1.0;

// catch the small flux factor case to prevent nans
if(fluxfac < minfluxfac)
Z = 3.*fluxfac;
else{
Real residual = 1.0;
int count = 0;
while(std::abs(residual)>maxresidual and count<maxcount){
residual = minerbo_residual(fluxfac, Z);
Real slope = minerbo_residual_derivative(fluxfac, Z);
Z -= residual/slope;
count++;
}
if(residual>maxresidual)
amrex::Error("Failed to converge on a solution.");
}

amrex::Print() << "fluxfac="<<fluxfac<<" Z=" << Z<<std::endl;
return Z;
}

namespace
{
AMREX_GPU_HOST_DEVICE void get_position_unit_cell(Real* r, const IntVect& nppc, int i_part)
Expand Down Expand Up @@ -69,6 +112,16 @@ namespace
*Usymmetric = 2. * (amrex::Random()-0.5);
}

// angular structure as determined by the Minerbo closure
// Z is a parameter determined by the flux factor
// mu is the cosine of the angle relative to the flux direction
// Coefficients set such that the expectation value is 1
AMREX_GPU_HOST_DEVICE void minerbo_closure(Real* result, const Real Z, const Real mu){
Real minfluxfac = 1e-3;
*result = std::exp(Z*mu);
if(Z/3.0 > minfluxfac)
dwillcox marked this conversation as resolved.
Show resolved Hide resolved
*result *= Z/std::sinh(Z);
dwillcox marked this conversation as resolved.
Show resolved Hide resolved
}
}

FlavoredNeutrinoContainer::
Expand Down Expand Up @@ -102,6 +155,27 @@ InitParticles(const TestParams* parms)

const Real scale_fac = dx[0]*dx[1]*dx[2]/nlocs_per_cell/ndirs_per_loc;

// get the Z parameters for the Minerbo closure if using simuation type 5
Real Ze, Za, Zx;
Real fluxfac_e, fluxfac_a, fluxfac_x;
if(parms->simulation_type==5){
fluxfac_e = std::sqrt(
parms->st5_fxnue*parms->st5_fxnue +
parms->st5_fynue*parms->st5_fynue +
parms->st5_fznue*parms->st5_fznue );
fluxfac_a = std::sqrt(
parms->st5_fxnua*parms->st5_fxnua +
parms->st5_fynua*parms->st5_fynua +
parms->st5_fznua*parms->st5_fznua );
fluxfac_x = std::sqrt(
parms->st5_fxnux*parms->st5_fxnux +
parms->st5_fynux*parms->st5_fynux +
parms->st5_fznux*parms->st5_fznux );
Ze = minerbo_Z(fluxfac_e);
Za = minerbo_Z(fluxfac_a);
Zx = minerbo_Z(fluxfac_x);
}

#ifdef _OPENMP
#pragma omp parallel
#endif
Expand Down Expand Up @@ -494,6 +568,82 @@ InitParticles(const TestParams* parms)
p.rdata(PIdx::Nbar) = ndensbar*scale_fac * (1. + 3.*parms->st4_fluxfacbar*costhetabar);
}

//====================//
// 5- Minerbo Closure //
//====================//
else if(parms->simulation_type==5){
AMREX_ASSERT(NUM_FLAVORS==3 or NUM_FLAVORS==2);

// set energy to 50 MeV
p.rdata(PIdx::pupt) = 50. * 1e6*CGSUnitsConst::eV;
p.rdata(PIdx::pupx) = u[0] * p.rdata(PIdx::pupt);
p.rdata(PIdx::pupy) = u[1] * p.rdata(PIdx::pupt);
p.rdata(PIdx::pupz) = u[2] * p.rdata(PIdx::pupt);

// get the cosine of the angle between the direction and each flavor's flux vector
Real mue = fluxfac_e>0 ? (parms->st5_fxnue*u[0] + parms->st5_fynue*u[1] + parms->st5_fznue*u[2])/fluxfac_e : 0;
Real mua = fluxfac_a>0 ? (parms->st5_fxnua*u[0] + parms->st5_fynua*u[1] + parms->st5_fznua*u[2])/fluxfac_a : 0;
Real mux = fluxfac_x>0 ? (parms->st5_fxnux*u[0] + parms->st5_fynux*u[1] + parms->st5_fznux*u[2])/fluxfac_x : 0;
dwillcox marked this conversation as resolved.
Show resolved Hide resolved

// get the number of each flavor in this particle.
// parms->st5_nnux contains the number density of mu+tau neutrinos+antineutrinos
// Nnux_thisparticle contains the number of EACH of mu/tau anti/neutrinos (hence the factor of 4)
Real angular_factor;
minerbo_closure(&angular_factor, Ze, mue);
Real Nnue_thisparticle = parms->st5_nnue*scale_fac * angular_factor;
minerbo_closure(&angular_factor, Za, mua);
Real Nnua_thisparticle = parms->st5_nnua*scale_fac * angular_factor;
minerbo_closure(&angular_factor, Zx, mux);
Real Nnux_thisparticle = parms->st5_nnux*scale_fac * angular_factor / 4.0;

// set total number of neutrinos the particle has as the sum of the flavors
p.rdata(PIdx::N ) = Nnue_thisparticle + Nnux_thisparticle;
p.rdata(PIdx::Nbar) = Nnua_thisparticle + Nnux_thisparticle;
#if NUM_FLAVORS==3
p.rdata(PIdx::N ) += Nnux_thisparticle;
p.rdata(PIdx::Nbar) += Nnux_thisparticle;
#endif

// set on-diagonals to have relative proportion of each flavor
p.rdata(PIdx::f00_Re) = Nnue_thisparticle / p.rdata(PIdx::N );
p.rdata(PIdx::f11_Re) = Nnux_thisparticle / p.rdata(PIdx::N );
p.rdata(PIdx::f00_Rebar) = Nnua_thisparticle / p.rdata(PIdx::Nbar);
p.rdata(PIdx::f11_Rebar) = Nnux_thisparticle / p.rdata(PIdx::Nbar);
#if NUM_FLAVORS==3
p.rdata(PIdx::f22_Re) = Nnux_thisparticle / p.rdata(PIdx::N );
p.rdata(PIdx::f22_Rebar) = Nnux_thisparticle / p.rdata(PIdx::Nbar);
#endif

// random perturbations to the off-diagonals
Real rand;
symmetric_uniform(&rand);
p.rdata(PIdx::f01_Re) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Re ) - p.rdata(PIdx::f11_Re ));
symmetric_uniform(&rand);
p.rdata(PIdx::f01_Im) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Re ) - p.rdata(PIdx::f11_Re ));
symmetric_uniform(&rand);
p.rdata(PIdx::f01_Rebar) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Rebar) - p.rdata(PIdx::f11_Rebar));
symmetric_uniform(&rand);
p.rdata(PIdx::f01_Imbar) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Rebar) - p.rdata(PIdx::f11_Rebar));
#if NUM_FLAVORS==3
symmetric_uniform(&rand);
p.rdata(PIdx::f02_Re) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Re ) - p.rdata(PIdx::f22_Re ));
symmetric_uniform(&rand);
p.rdata(PIdx::f02_Im) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Re ) - p.rdata(PIdx::f22_Re ));
symmetric_uniform(&rand);
p.rdata(PIdx::f12_Re) = parms->st5_amplitude*rand * (p.rdata(PIdx::f11_Re ) - p.rdata(PIdx::f22_Re ));
symmetric_uniform(&rand);
p.rdata(PIdx::f12_Im) = parms->st5_amplitude*rand * (p.rdata(PIdx::f11_Re ) - p.rdata(PIdx::f22_Re ));
symmetric_uniform(&rand);
p.rdata(PIdx::f02_Rebar) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Rebar) - p.rdata(PIdx::f22_Rebar));
symmetric_uniform(&rand);
p.rdata(PIdx::f02_Imbar) = parms->st5_amplitude*rand * (p.rdata(PIdx::f00_Rebar) - p.rdata(PIdx::f22_Rebar));
symmetric_uniform(&rand);
p.rdata(PIdx::f12_Rebar) = parms->st5_amplitude*rand * (p.rdata(PIdx::f11_Rebar) - p.rdata(PIdx::f22_Rebar));
symmetric_uniform(&rand);
p.rdata(PIdx::f12_Imbar) = parms->st5_amplitude*rand * (p.rdata(PIdx::f11_Rebar) - p.rdata(PIdx::f22_Rebar));
#endif
}

else{
amrex::Error("Invalid simulation type");
}
Expand Down
23 changes: 23 additions & 0 deletions Source/Parameters.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,13 @@ struct TestParams : public amrex::Gpu::Managed
Real st4_ndensbar, st4_thetabar, st4_phibar, st4_fluxfacbar;
Real st4_amplitude;

// simulation_type==5
Real st5_nnue , st5_nnua , st5_nnux ;
Real st5_fxnue, st5_fxnua, st5_fxnux;
Real st5_fynue, st5_fynua, st5_fynux;
Real st5_fznue, st5_fznua, st5_fznux;
Real st5_amplitude;

void Initialize(){
ParmParse pp;
pp.get("simulation_type", simulation_type);
Expand Down Expand Up @@ -105,6 +112,22 @@ struct TestParams : public amrex::Gpu::Managed
pp.get("st4_fluxfacbar", st4_fluxfacbar);
pp.get("st4_amplitude", st4_amplitude);
}

if(simulation_type==5){
pp.get("st5_nnue",st5_nnue);
pp.get("st5_nnua",st5_nnua);
pp.get("st5_nnux",st5_nnux);
pp.get("st5_fxnue",st5_fxnue);
pp.get("st5_fxnua",st5_fxnua);
pp.get("st5_fxnux",st5_fxnux);
pp.get("st5_fynua",st5_fynua);
pp.get("st5_fynux",st5_fynux);
pp.get("st5_fynue",st5_fynue);
pp.get("st5_fznue",st5_fznue);
pp.get("st5_fznua",st5_fznua);
pp.get("st5_fznux",st5_fznux);
pp.get("st5_amplitude",st5_amplitude);
}
}
};

Expand Down