Skip to content

Commit

Permalink
added 1D and 2D oblique shock profiles
Browse files Browse the repository at this point in the history
  • Loading branch information
Sophie Aerdker committed Oct 4, 2023
1 parent d091d2c commit 7daa74e
Show file tree
Hide file tree
Showing 2 changed files with 303 additions and 6 deletions.
98 changes: 92 additions & 6 deletions include/crpropa/advectionField/AdvectionField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,96 @@ class SphericalAdvectionField: public AdvectionField {
std::string getDescription() const;
};

/**
@class 1DAdvectionShock
@brief Advection field in x-direction with shock at x = 0 and width x_sh approximated by tanh()
with variable compression ratio r_comp = v_up/v_down
*/
class OneDimensionalAdvectionShock: public AdvectionField {
double r_comp; //compression ratio of shock
double v_up; //upstream velocity
double x_sh; //shock width
public:/** Constructor
@param r_comp //compression ratio of shock
@param v_up //upstream velocity
@param x_sh //shock width
*/
OneDimensionalAdvectionShock(double r_comp, double v_up, double x_sh);
Vector3d getField(const Vector3d &position) const;
double getDivergence(const Vector3d &position) const;

void setComp(double r_comp);
void setVup(double v_up);
void setShockwidth(double x_sh);


std::string getDescription() const;
};

/**
@class 1DSphericalShock
@brief Advection field in r-direction with shock at r_sh and width l_sh approximated by tanh()
with variable compression ratio r_comp = v_up/v_down
*/
class OneDimensionalSphericalShock: public AdvectionField {
double r_comp; //compression ratio of shock
double v_up; //upstream velocity
double l_sh; //shock width
double r_sh; //shock radius
bool cool_ups; //flag for upstream cooling
public:/** Constructor
@param r_comp //compression ratio of shock
@param v_up //upstream velocity
@param l_sh //shock width
@param r_sh //shock radius
@param cool_ups //flag for upstream cooling
*/
OneDimensionalSphericalShock(double r_sh, double v_up, double q, double l_sh, bool cool_ups);
Vector3d getField(const Vector3d &position) const;
double getDivergence(const Vector3d &position) const;

void setComp(double q);
void setVup(double v_up);
void setShockwidth(double l_sh);
void setShockRadius(double r_sh);
void setCooling(bool cool_ups);

std::string getDescription() const;
};

/**
@class ObliqueAdvectionShock
@brief Advection field in x-y-direction with shock at x = 0 and width x_sh approximated by tanh()
with variable compression ratio r_comp = vx_up/vx_down. The y component vy is not shocked
and remains constant.
*/
class ObliqueAdvectionShock: public AdvectionField {
double r_comp; //compression ratio of shock
double vx_up; //upstream velocity x-component
double vy; //constant velocity y-component
double x_sh; //shock width

public:/** Constructor
@param r_comp //compression ratio of shock
@param vx_up //upstream velocity x-component
@param vy //constant velocity y-component
@param x_sh //shock width
*/
ObliqueAdvectionShock(double r_comp, double vx_up, double vy, double x_sh);
Vector3d getField(const Vector3d &position) const;
double getDivergence(const Vector3d &position) const;

void setComp(double r_comp);
void setVup(double vx_up);
void setVy(double vy);
void setShockwidth(double x_sh);


std::string getDescription() const;
};

/**
/**
@class SphericalAdvectionShock
Expand All @@ -150,6 +240,8 @@ class SphericalAdvectionShock: public AdvectionField {
@param r_0 Position of the shock
@param v_0 Constant velocity (r<<r_o)
@param lambda Transition width / width of the shock
@param r_rot Normalization radius for rotation speed
@param v_phi Rotation speed at r_rot
*/
SphericalAdvectionShock(const Vector3d origin, double r_0, double v_0, double lambda);

Expand All @@ -163,13 +255,7 @@ class SphericalAdvectionShock: public AdvectionField {
void setR0(double r);
void setV0(double v);
void setLambda(double l);
/**
* @param r Normalization radius for rotation speed
*/
void setRRot(double r);
/**
* @param vPhi Rotation speed at r_rot
*/
void setAzimuthalSpeed(double vPhi);

Vector3d getOrigin() const;
Expand Down
211 changes: 211 additions & 0 deletions src/advectionField/AdvectionField.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,217 @@ std::string SphericalAdvectionField::getDescription() const {
return s.str();
}

//----------------------------------------------------------------

OneDimensionalAdvectionShock::OneDimensionalAdvectionShock(double r_comp, double v_up, double x_sh){

setComp(r_comp);
setVup(v_up);
setShockwidth(x_sh);
}

Vector3d OneDimensionalAdvectionShock::getField(const Vector3d &position) const {

double x = position.x;
double v_down = v_up/r_comp;

double a = (v_up + v_down)*0.5;
double b = (v_up - v_down)*0.5;

Vector3d v(0.);
v.x = a - b*tanh(x/x_sh);

return v;
}

double OneDimensionalAdvectionShock::getDivergence(const Vector3d &position) const {

double x = position.x;
double v_down = v_up/r_comp;

double a = (v_up + v_down)*0.5;
double b = (v_up - v_down)*0.5;

return -b/x_sh*(1 - tanh(x/x_sh)*tanh(x/x_sh));
}

void OneDimensionalAdvectionShock::setComp(double r) {
r_comp = r;
return;
}

void OneDimensionalAdvectionShock::setVup(double v) {
v_up = v;
return;
}
void OneDimensionalAdvectionShock::setShockwidth(double w) {
x_sh = w;
return;
}

std::string OneDimensionalAdvectionShock::getDescription() const {
std::stringstream s;
s << "Shock width: " << x_sh / km << " km, ";
s << "Vup: " << v_up / km * sec << " km/s, ";
s << "Comp: " << r_comp;

return s.str();
}

//----------------------------------------------------------------

OneDimensionalSphericalShock::OneDimensionalSphericalShock(double r_sh, double v_up, double q, double l_sh, bool cool_ups){

setComp(q);
setVup(v_up);
setShockwidth(l_sh);
setShockRadius(r_sh);
setCooling(cool_ups);
}

Vector3d OneDimensionalSphericalShock::getField(const Vector3d &position) const {

double r = position.getR();
Vector3d e_r = position.getUnitVector();

double v_down = v_up/r_comp;
double a = (v_up + v_down)*0.5;
double b = (v_up - v_down)*0.5;

double v;
if (cool_ups == true){

if (r <= r_sh)
v = a - b*tanh((r-r_sh)/l_sh);
else
v = (a - b*tanh((r-r_sh)/l_sh)) * (r_sh/r)*(r_sh/r);

}
else
v = (a - b*tanh((r-r_sh)/l_sh)) * (r_sh/r)*(r_sh/r);

return v * e_r;
}

double OneDimensionalSphericalShock::getDivergence(const Vector3d &position) const {

double r = position.getR();

double v_down = v_up/r_comp;

double a = (v_up + v_down)*0.5;
double b = (v_up - v_down)*0.5;

if (cool_ups == true){
if (r <= r_sh)
return 2*a/r - 2*b/r * tanh((r-r_sh)/l_sh) - b/l_sh * (1 - tanh((r-r_sh)/l_sh)*tanh((r-r_sh)/l_sh));
else
return -(r_sh/r)*(r_sh/r) * b/l_sh * (1 - tanh((r-r_sh)/l_sh)*tanh((r-r_sh)/l_sh));
}
else
return -(r_sh/r)*(r_sh/r) * b/l_sh * (1 - tanh((r-r_sh)/l_sh)*tanh((r-r_sh)/l_sh));

}

void OneDimensionalSphericalShock::setComp(double r) {
r_comp = r;
return;
}

void OneDimensionalSphericalShock::setVup(double v) {
v_up = v;
return;
}
void OneDimensionalSphericalShock::setShockwidth(double w) {
l_sh = w;
return;
}

void OneDimensionalSphericalShock::setShockRadius(double r) {
r_sh = r;
return;
}

void OneDimensionalSphericalShock::setCooling(bool c) {
cool_ups = c;
return;
}

std::string OneDimensionalSphericalShock::getDescription() const {
std::stringstream s;
s << "Shock width: " << l_sh / km << " km, ";
s << "Shock radius: " << r_sh / km << " km, ";
s << "Vup: " << v_up / km * sec << " km/s, ";
s << "Comp: " << r_comp;

return s.str();
}

//----------------------------------------------------------------

ObliqueAdvectionShock::ObliqueAdvectionShock(double r_comp, double vx_up, double vy, double x_sh){

setComp(r_comp);
setVup(vx_up);
setVy(vy);
setShockwidth(x_sh);
}

Vector3d ObliqueAdvectionShock::getField(const Vector3d &position) const {

double x = position.x;
double vx_down = vx_up/r_comp;

double a = (vx_up + vx_down)*0.5;
double b = (vx_up - vx_down)*0.5;

Vector3d v(0.);
v.x = a - b*tanh(x/x_sh);
v.y = vy;

return v;
}

double ObliqueAdvectionShock::getDivergence(const Vector3d &position) const {

double x = position.x;
double vx_down = vx_up/r_comp;
// vy = const

double a = (vx_up + vx_down)*0.5;
double b = (vx_up - vx_down)*0.5;

return -b/x_sh*(1 - tanh(x/x_sh)*tanh(x/x_sh));
}

void ObliqueAdvectionShock::setComp(double r) {
r_comp = r;
return;
}

void ObliqueAdvectionShock::setVup(double v) {
vx_up = v;
return;
}

void ObliqueAdvectionShock::setVy(double v) {
vy = v;
return;
}
void ObliqueAdvectionShock::setShockwidth(double w) {
x_sh = w;
return;
}

std::string ObliqueAdvectionShock::getDescription() const {
std::stringstream s;
s << "Shock width: " << x_sh / km << " km, ";
s << "Vx_up: " << vx_up / km * sec << " km/s, ";
s << "Vy: " << vy / km * sec << " km/s, ";
s << "Comp: " << r_comp;

return s.str();
}

//-----------------------------------------------------------------

Expand Down

0 comments on commit 7daa74e

Please sign in to comment.