Homec4science
Diffusion Gear (master) Browse

Gear/src/outerpotential.c
b87f34a0ce7emaster

Last Change
Raw File

outerpotential.c
View Options

#include <dirent.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include "allvars.h"
#include "proto.h"
/**************************************
*
* * * * driving functions * * *
*
***************************************/
#ifdef OUTERPOTENTIAL
void init_outer_potential()
{
#ifdef NFW
init_outer_potential_nfw();
#endif
#ifdef PLUMMER
init_outer_potential_plummer();
#endif
#ifdef HOMOSPHERE
init_outer_potential_homosphere();
#endif
#ifdef MIYAMOTONAGAI
init_outer_potential_miyamotonagai();
#endif
#ifdef PISOTHERM
init_outer_potential_pisotherm();
#endif
#ifdef CORIOLIS
init_outer_potential_coriolis();
#endif
#ifdef EVOLVE_POT
init_outer_potential_evolve();
#endif
}
void outer_forces()
{
int i,j;
#ifdef EVOLVE_POT
interp_evolve_vals();
#endif
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
FLOAT pos[3] = {}, vel[3] = {}, acc[3] = {};
#ifdef PERIODICOUTER //The position is (the internal position - the origin of the tracer particle)*rotated to the galactic axis + tracer particle
double posquat[4] = {P[i].Pos[0] - All.TracePeriodPos[0],P[i].Pos[1] - All.TracePeriodPos[1],P[i].Pos[2] - All.TracePeriodPos[2],0};
double velquat[4] = {P[i].Vel[0],P[i].Vel[1],P[i].Vel[2],0};
quatrot(posquat,posquat,All.TraceQuat);
quatrot(velquat,velquat,All.TraceQuat);
for(j = 0; j<3; j++) //Last entry of posbase should be 0, so can ignore
{
pos[j] = posquat[j] + All.TracePos[j] ; //Having it have a position makes it easier to calculate initial velocity
vel[j] = velquat[j] + All.TraceVel[j];
}
#else //Position is just the position
for( j = 0; j<3; j++)
{
pos[j] = P[i].Pos[j];
vel[j] = P[i].Vel[j];
}
#endif
forces_potentials(pos,vel,acc); //Calculate acceleration
#ifdef PERIODICOUTER //Remove acceleration of the tracer
FLOAT rotpos[3];
double accquat[4] = {};
for( j = 0; j<3; j++)
{
acc[j] -= All.TraceAcc[j];
accquat[j] = acc[j];
}
quatrot(accquat,accquat,All.TraceQuatConj); //Rotate back
for( j=0;j<3;j++)
{
acc[j] = accquat[j];
}
#endif
for( j = 0; j<3; j++)
{
P[i].GravAccel[j] += acc[j]; //Add final acceleration to the gravaccel
}
}
}
}
void forces_potentials(FLOAT pos[3],FLOAT vel[3],FLOAT acc[3])
{
#ifdef NFW
outer_forces_nfw(pos,acc);
#endif
#ifdef PLUMMER
outer_forces_plummer(pos,acc);
#endif
#ifdef HOMOSPHERE
outer_forces_homosphere(pos,acc);
#endif
#ifdef PISOTHERM
outer_forces_pisotherm(pos,acc);
#endif
#ifdef MIYAMOTONAGAI
outer_forces_miyamotonagai(pos,acc);
#endif
#ifdef CORIOLIS
outer_forces_coriolis(pos,vel,acc);
#endif
#ifdef EVOLVE_POT
outer_forces_evolve(pos,acc);
#endif
}
void outer_potential()
{
int i,j;
#ifdef EVOLVE_POT
interp_evolve_vals();
#endif
for (i= 0; i < NumPart; i++)
{
FLOAT pos[3];
FLOAT Phi[1] = {};
#ifdef PERIODICOUTER //Add the tracer particle position to get the real radius
double posquat[4] = {P[i].Pos[0] - All.TracePeriodPos[0],P[i].Pos[1] - All.TracePeriodPos[1],P[i].Pos[2] - All.TracePeriodPos[2],0};
quatrot(posquat,posquat,All.TraceQuat);
for(j=0 ; j<3; j++)
{
pos[j] = posquat[j] + All.TracePos[j];
}
#else
for(j=0 ; j<3; j++)
{
pos[j] = P[i].Pos[j];
}
#endif
potential_components(pos,Phi);
/* !!! do not forget do multiply by 2 (Springel choice ?) !!! */
P[i].Potential = P[i].Potential + 2*Phi[0];
}
}
void potential_components(FLOAT pos[3], FLOAT Phi[1])
{
#ifdef NFW
outer_potential_nfw(pos,Phi);
#endif
#ifdef PLUMMER
outer_potential_plummer(pos,Phi);
#endif
#ifdef HOMOSPHERE
outer_potential_homosphere(pos,Phi);
#endif
#ifdef MIYAMOTONAGAI
outer_potential_miyamotonagai(pos,Phi);
#endif
#ifdef PISOTHERM
outer_potential_pisotherm(pos,Phi);
#endif
#ifdef CORIOLIS
outer_potential_coriolis(pos,Phi);
#endif
#ifdef EVOLVE_POT
outer_potential_evolve(pos,Phi);
#endif
}
#ifdef NFW
/****************************
*
* * * * NFW * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_nfw()
{
double rhoc;
All.HaloMass = All.HaloMass / All.UnitMass_in_g;
All.HaloMass = All.HaloMass*SOLAR_MASS; /* from Msolar in user units */
All.NFWPotentialCte = (1.-All.GasMassFraction) * All.G*All.HaloMass/(log(1+All.HaloConcentration) - All.HaloConcentration/(1+All.HaloConcentration));
/* here, convert HubbleParam in 1/(unit_time) */
rhoc = pow(All.HubbleParam*HUBBLE*All.UnitTime_in_s,2)*3/(8*PI*All.G);
All.Rs = pow(( All.HaloMass/(4/3.*PI*200*rhoc) ),(1./3.)) / All.HaloConcentration;
}
/*! compute forces du to the outer potential
*
*/
void outer_forces_nfw(FLOAT pos[3], FLOAT acc[3])
{
FLOAT dPhi, r;
int j;
r = sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
dPhi = -All.NFWPotentialCte * ( log(1+r/All.Rs) - r/(All.Rs+r) ) / (r*r);
for(j = 0; j < 3; j++)
{
acc[j] += dPhi * pos[j]/r;
}
}
/*! compute outer potential
*
*/
void outer_potential_nfw(FLOAT pos[3], FLOAT Phi[1])
{
FLOAT r;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
Phi[0] += - All.NFWPotentialCte * (log(1+r/All.Rs) / r);
}
#endif
#ifdef PLUMMER
/****************************
*
* * * * PLUMMER * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_plummer()
{
All.PlummerPotentialCte = All.G*All.PlummerMass /All.UnitMass_in_g *SOLAR_MASS; /* from Msolar in user units */
}
/*! compute forces due to the outer potential
*
*/
void outer_forces_plummer(FLOAT pos[3], FLOAT acc[3])
{
FLOAT r,dPhi;
int j;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
dPhi = - All.PlummerPotentialCte / pow(r*r+All.PlummerSoftenning*All.PlummerSoftenning,1.5)*r;
for(j = 0; j<3; j++)
{
acc[j] += dPhi * pos[j]/r;
}
}
/*! compute outer potential
*
*/
void outer_potential_plummer(FLOAT pos[3], FLOAT Phi[1])
{
FLOAT r;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
Phi[0] += - All.PlummerPotentialCte / pow(r*r+All.PlummerSoftenning*All.PlummerSoftenning,0.5);
}
#endif
#ifdef HOMOSPHERE
/***************************************
*
* * * * HOMOGENEOUS SPHERE * * *
*
****************************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_homosphere()
{
//All.HomosphereMass = All.G*All.HomosphereMass /All.UnitMass_in_g *SOLAR_MASS; /* from Msolar in user units */
All.HomospherePotentialCte = All.G*All.HomosphereMass;
}
/*! compute forces due to the outer potential
*
*/
void outer_forces_homosphere(FLOAT pos[3], FLOAT acc[3])
{
FLOAT r,dPhi;
int j;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
if (r <= All.HomosphereMaxRadius)
dPhi = -All.HomospherePotentialCte / pow(All.HomosphereMaxRadius,3) * r;
else
dPhi = -All.HomospherePotentialCte / pow(r,2);
for(j = 0; j<3; j++)
{
acc[j] += dPhi * pos[j]/r;
}
}
/*! compute outer potential
*
*/
void outer_potential_homosphere(FLOAT pos[3], FLOAT Phi[1])
{
FLOAT r;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
if (r <= All.HomosphereMaxRadius)
Phi[0] += - All.HomospherePotentialCte / pow(All.HomosphereMaxRadius,3) * pow(r,2);
else
Phi[0] += - All.HomospherePotentialCte / r;
}
#endif
#ifdef MIYAMOTONAGAI
/****************************
*
* * * * MIYAMOTONAGAI * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_miyamotonagai()
{
All.MiyamotoNagaiPotentialCte = All.G*All.MiyamotoNagaiMass /All.UnitMass_in_g *SOLAR_MASS; /* from Msolar in user units */
}
/*! compute forces due to the outer potential
*
*/
void outer_forces_miyamotonagai(FLOAT pos[3], FLOAT acc[3])
{
FLOAT r,dPhi;
int j;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
dPhi = - All.MiyamotoNagaiPotentialCte / pow(r*r+All.MiyamotoNagaiHr*All.MiyamotoNagaiHr,1.5)*r;
for(j = 0; j<3; j++)
{
acc[j] += dPhi * pos[j]/r;
}
}
/*! compute outer potential
*
*/
void outer_potential_miyamotonagai(FLOAT pos[3], FLOAT Phi[1])
{
FLOAT r;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
Phi[0] += - All.MiyamotoNagaiPotentialCte / pow(r*r+All.MiyamotoNagaiHr*All.MiyamotoNagaiHr,0.5);
}
#endif
#ifdef PISOTHERM
/**********************************************
*
* * * * PSEUDO ISOTHERM POTENTIAL * * *
*
***********************************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_pisotherm()
{
All.PisothermPotentialCte = 4*PI*All.Rho0*pow(All.Rc,3) *All.G *(1.-All.GasMassFraction);
All.Potentialw = gsl_integration_workspace_alloc (1000);
All.PotentialF.function = &potential_f;
/* potential normalisation */
/*All.PotentialInf = All.PisothermPotentialCte/All.Rc; viriel Pfen. */
All.PotentialInf = 3.3263051716357208; /* in order to have 2Uint = Epot */
}
/*! potential function
*
*/
double potential_f (double r, void * params)
{
double f = All.PisothermPotentialCte*(r/All.Rc-atan(r/All.Rc))/(r*r);
return f;
}
/*! potential function
*
*/
double get_potential(double r)
{
double result, error;
gsl_integration_qags (&All.PotentialF, 0, r, 0, 1e-5, 1000, All.Potentialw, &result, &error);
return result;
}
/*! compute forces du to the outer potential
*
*/
void outer_forces_pisotherm(FLOAT pos[3], FLOAT acc[3])
{
FLOAT r,dPhi;
int j;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
dPhi = - All.PisothermPotentialCte * ( r/All.Rc - atan(r/All.Rc) ) / (r*r);
for(j = 0; j < 3; j++)
{
acc[j] += dPhi * pos[j]/r;
}
}
/*! compute outer potential
*
*/
void outer_potential_pisotherm(FLOAT pos[3], FLOAT Phi[1])
{
FLOAT r;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
Phi[0] += get_potential(r)-All.PotentialInf;
}
#endif
#ifdef CORIOLIS
/****************************
*
* * * * CORIOLIS * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_coriolis()
{
set_outer_potential_coriolis();
}
void set_outer_potential_coriolis()
{
double amp;
double CoriolisTime = 1000;
amp = All.Time/CoriolisTime;
if (amp>1)
amp = 1;
All.CoriolisOmegaX = All.CoriolisOmegaX0*amp;
All.CoriolisOmegaY = All.CoriolisOmegaY0*amp;
All.CoriolisOmegaZ = All.CoriolisOmegaZ0*amp;
printf("coriolis CoriolisOmegaZ = %g\n",All.CoriolisOmegaZ);
}
/*! compute forces due to the outer potential
*
*/
void outer_forces_coriolis(FLOAT pos[3], FLOAT acc[3])
{
acc[0] += pow(All.CoriolisOmegaZ,2)*pos[0] - All.CoriolisOmegaZ*P[i].Vel[1];
acc[1] += pow(All.CoriolisOmegaZ,2)*pos[1] + All.CoriolisOmegaZ*P[i].Vel[0];
}
/*! compute outer potential
*
*/
void outer_potential_coriolis(FLOAT pos[3], FLOAT Phi[1])
{
FLOAT r;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
Phi[0] += -0.5*pow(All.CoriolisOmegaZ*r,2);
}
#endif
#ifdef EVOLVE_POT
/*NOTE THE FOLLOWING ASSUMES A LINUX FILE STRUCTURE*/
/*Since we can't broadcast splines we have two choices. 1) Make every CPU calculate the spline 2) At every timestep, broadcast information evaluated from the spline. We follow option 1 here.*/
/*Error checking is pretty terrible, i.e. nearly non-existent*/
void init_outer_potential_evolve()
{
struct dirent **potnames;
int npots, i, j,k,nlines,nheader;
int LINELENGTH = 8192;
char pfilename[MAXLEN_FILENAME];
char line[LINELENGTH];
char pottype[POTNAMELENGTH]; //pottype decides on the type of profile. 0 = NFW, 1 = Plummer.
FILE *potential;
FILE *potentialline;
double entry;
int nvars;
int coeffs = 12;
int order = 4;
/*Read folder. Count total number of entries. Create array with that many.*/
npots = scandir(All.PotDirect, &potnames,0,0);
/* ./ and ../ aren't potential files, so can free them automatically and remove them from the count, also if the files are hidden well its their own fdamn fault*/
free(potnames[0]);
free(potnames[1]);
npots -= 2;
if(npots > MAXPOTFILE)
{
if(ThisTask == 0)
printf("\nNumber of potentials %d > MAXPOTFILE %d\n",npots,MAXPOTFILE);\
endrun(5000);
}
/*For each potential file, need to interpolate all the variables used throughout, e.g. Mvir, Rvir, Rs for NFW.*/
for(i=0;i<npots;i++)
{
sprintf(pfilename,"%s/%s",All.PotDirect,potnames[i+2]->d_name);
/*Open the potential file*/
potential = fopen(pfilename,"r");
free(potnames[i+2]);
nlines = 0;
nheader = 0;
while(fgets(line,LINELENGTH,potential) != NULL)
{
/*Allow headers to exist but only at the start*/
if(nlines == 0 && &line[0] == "#")
nheader++;
else
nlines++;
}
rewind(potential); /*Have the number of lines, so move back to start of file*/
nlines--; /*First non header line is the type so don't want to count it */
for(j=0;j<nheader;j++)
fgets(line,LINELENGTH,potential); /*Grab each header line so we can skip it...*/
fgets(line,LINELENGTH,potential); /*First line*/
sscanf(line,"%s",pottype);
if(ThisTask == 0)
printf("%s\n",pottype);
All.PotType[i] = pottype;
///*How many variables are we expecting and setting potential type, don't change these just add more*/
//if(!strcasecmp(pottype,"NFW")) All.PotType[i] = 0;
//else if(!strcasecmp(pottype,"Plummer")) All.PotType[i] = 1;
nvars = NumVars(All.PotType[i]);
double vars[nvars+1][nlines];
/*Grab each lines variables and throw them into the vars array*/
/*Move to a file to use fscanf, since various files might have different lines*/
for(j=0;j<nlines;j++)
{
fgets(line,LINELENGTH,potential);
potentialline = fmemopen(line,strlen(line),"r");
for(k=0;k<nvars+1;k++)
fscanf(potentialline,"%lg",&vars[k][j]);
fclose(potentialline);
}
fclose(potential);
adapt_evolve_variables(All.PotType[i],*vars,nlines);
/*Interpolate each of these with respect to time which is the first column*/
for(j=0;j<nvars;j++)
{
All.PotSplines[i][j] = gsl_spline_alloc (gsl_interp_cspline, nlines);
gsl_spline_init(All.PotSplines[i][j],vars[0],vars[j+1],nlines);
}
}
free(potnames);
for(i=npots;i<MAXPOTFILE;i++)
{
All.PotType[i] = -1; //Casting int to char but it's fine...
}
interp_evolve_vals();
}
int NumVars(char pottype[POTNAMELENGTH])
{
if(pottype == "NFW") //NFW
return 3;
else if (pottype=="Plummer") //Plummer
return 2;
}
void adapt_evolve_variables(char pottype[POTNAMELENGTH], double *vars, int nlines)
{
int i,nvars;
if(pottype == "NFW")
for(i=0;i<nlines;i++)
{
nvars = 4;
*((vars+i)) *= SEC_PER_MEGAYEAR*1e3/All.UnitTime_in_s*All.HubbleParam; //Time in Gyr
*((vars+nlines+i)) *= SOLAR_MASS/All.UnitMass_in_g*All.HubbleParam; //Mvir in Msun
*((vars+2*nlines+i)) *= CM_PER_MPC/1e3/All.UnitLength_in_cm*All.HubbleParam; //Rvir in kpc
*((vars+3*nlines+i)) *= CM_PER_MPC/1e3/All.UnitLength_in_cm*All.HubbleParam; //Rs in kpc
}
else if(pottype == "Plummer")
for(i=0;i<nlines;i++)
{
nvars = 3;
*((vars+i)) *= 1000./All.UnitTime_in_Megayears*All.HubbleParam; //Time in Gyr
*((vars+nlines+i)) *= SOLAR_MASS/All.UnitMass_in_g*All.HubbleParam; //Mass in Msun
*((vars+2*nlines+i)) *= 30.857e20/All.UnitLength_in_cm*All.HubbleParam; //Softening in kpc
}
}
void interp_evolve_vals()
{
gsl_interp_accel *spline_acc = gsl_interp_accel_alloc();
int i,j,nvars;
double val[MAXNVARS];
double work[5]; //Values to use for temp calculations need to be commented what they do
for(i=0;i<MAXPOTFILE;i++)
{
if(All.PotType[i] == -1)
break;
gsl_interp_accel_reset(spline_acc);
nvars = NumVars(All.PotType[i]);
for(j=0;j<nvars;j++)
val[j] = gsl_spline_eval(All.PotSplines[i][j],All.Time,spline_acc);
if(All.PotType[i] == "NFW")
{
work[0] = val[1]/val[2]; //Concentration
All.PotVals[i][0] = All.G*val[0]/(log(1+work[0]) - work[0]/(1+work[0]));//All.NFWPotentialCte Equivalent
All.PotVals[i][1] = val[2];
}
else if(All.PotType[i] == "Plummer")
{
All.PotVals[i][0] = All.G*val[0];
All.PotVals[i][1] = val[1];
}
}
gsl_interp_accel_free(spline_acc);
}
void outer_potential_evolve(FLOAT pos[3], FLOAT Phi[1])
{
FLOAT r;
int i;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
for(i=0;i<MAXPOTFILE;i++)
{
if(All.PotType[i] == -1)
break;
if(All.PotType[i] == "NFW")
{
Phi[0] += -All.PotVals[i][0] * (log(1+r/All.PotVals[i][1])/r);
}
else if(All.PotType[i] == "Plummer")
{
Phi[0] += -All.PotVals[i][0] / pow(r*r + All.PotVals[i][1]*All.PotVals[i][1],0.5);
}
}
}
void outer_forces_evolve(FLOAT pos[3], FLOAT acc[3])
{
FLOAT r,dPhi;
int i,j;
r = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
for(i=0;i<MAXPOTFILE;i++)
{
if(All.PotType[i] == -1)
break;
if(All.PotType[i] == "NFW")
{
dPhi = - All.PotVals[i][0] * ( log(1+r/All.PotVals[i][1]) - r/(All.PotVals[i][1]+r) ) / (r*r);
}
else if(All.PotType[i] == "Plummer")
{
dPhi = - All.PotVals[i][0] * ( r/All.PotVals[i][1] - atan(r/All.PotVals[i][1]) ) / (r*r);
}
for(j = 0; j < 3; j++)
{
acc[j] += dPhi * pos[j]/r;
}
}
}
#endif
#endif
c4science · Help