Skip to content

Commit

Permalink
Merge pull request #441 from sophieaerdker/ShockAdvFields
Browse files Browse the repository at this point in the history 
New advection fields
  • Loading branch information
@JulienDoerner
JulienDoerner authored Oct 16, 2023
2 parents cab3867 + 6cab22e commit 026c6d4
Show file tree
Hide file tree
Showing 2 changed files with 359 additions and 6 deletions.
116 changes: 110 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,110 @@ class SphericalAdvectionField: public AdvectionField {
std::string getDescription() const;
};

/**
@class OneDimensionalCartesianShock
@brief Advection field in x-direction with shock at x = 0 and width lShock approximated by tanh()
with variable compression ratio vUp/vDown
*/
class OneDimensionalCartesianShock: public AdvectionField {
double compressionRatio; //compression ratio of shock
double vUp; //upstream velocity
double lShock; //shock width
public:
/** Constructor
@param compressionRatio //compression ratio of shock
@param vUp //upstream velocity
@param lShock //shock width
*/
OneDimensionalCartesianShock(double compressionRatio, double vUp, double lShock);
Vector3d getField(const Vector3d &position) const;
double getDivergence(const Vector3d &position) const;

void setComp(double compressionRatio);
void setVup(double vUp);
void setShockwidth(double lShock);

double getComp() const;
double getVup() const;
double getShockwidth() const;

std::string getDescription() const;
};

/**
@class OneDimensionalSphericalShock
@brief Advection field in x-direction with shock at rShock and width lShock approximated by tanh()
with variable compression ratio ratio vUp/vDown
*/
class OneDimensionalSphericalShock: public AdvectionField {
double compressionRatio; //compression ratio of shock
double vUp; //upstream velocity
double lShock; //shock width
double rShock; //shock radius
bool coolUpstream; //flag for upstream cooling
public:
/** Constructor
@param compressionRatio //compression ratio of shock
@param vUp //upstream velocity
@param lShock //shock width
@param rShock //shock radius
@param coolUpstream //flag for upstream cooling
*/
OneDimensionalSphericalShock(double rShock, double vUp, double compressionRatio, double lShock, bool coolUpstream);
Vector3d getField(const Vector3d &position) const;
double getDivergence(const Vector3d &position) const;

void setComp(double compressionRatio);
void setVup(double vUp);
void setShockwidth(double lShock);
void setShockRadius(double rShock);
void setCooling(bool coolUpstream);

double getComp() const;
double getVup() const;
double getShockwidth() const;
double getShockRadius() const;
bool getCooling() const;

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 compressionRatio; //compression ratio of shock
double vXUp; //upstream velocity x-component
double vY; //constant velocity y-component
double lShock; //shock width

public:
/** Constructor
@param compressionRatio //compression ratio of shock
@param vXUp //upstream velocity x-component
@param vY //constant velocity y-component
@param lShock //shock width
*/
ObliqueAdvectionShock(double compressionRatio, double vXUp, double vY, double lShock);
Vector3d getField(const Vector3d &position) const;
double getDivergence(const Vector3d &position) const;

void setComp(double compressionRatio);
void setVup(double vXUp);
void setVy(double vY);
void setShockwidth(double lShock);

double getComp() const;
double getVup() const;
double getVy() const;
double getShockwidth() const;

std::string getDescription() const;
};

/**
@class SphericalAdvectionShock
Expand Down Expand Up @@ -163,20 +267,20 @@ 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;
double getR0() const;
double getV0() const;
double getLambda() const;
/**
* @param r Normalization radius for rotation speed
*/
double getRRot() const;
/**
* @param vPhi Rotation speed at r_rot
*/
double getAzimuthalSpeed() const;

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

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

OneDimensionalCartesianShock::OneDimensionalCartesianShock(double compressionRatio, double vUp, double lShock){
setComp(compressionRatio);
setVup(vUp);
setShockwidth(lShock);
}

Vector3d OneDimensionalCartesianShock::getField(const Vector3d &position) const {
double x = position.x;
double vDown = vUp / compressionRatio;

double a = (vUp + vDown) * 0.5;
double b = (vUp - vDown) * 0.5;

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

}

double OneDimensionalCartesianShock::getDivergence(const Vector3d &position) const {
double x = position.x;
double vDown = vUp / compressionRatio;

double a = (vUp + vDown) * 0.5;
double b = (vUp - vDown) * 0.5;
return -b / lShock * (1 - tanh(x / lShock) * tanh(x / lShock));
}

void OneDimensionalCartesianShock::setComp(double r) {
compressionRatio = r;
return;
}

void OneDimensionalCartesianShock::setVup(double v) {
vUp = v;
return;
}
void OneDimensionalCartesianShock::setShockwidth(double w) {
lShock = w;
return;
}

double OneDimensionalCartesianShock::getComp() const {
return compressionRatio;
}

double OneDimensionalCartesianShock::getVup() const {
return vUp;
}

double OneDimensionalCartesianShock::getShockwidth() const {
return lShock;
}

std::string OneDimensionalCartesianShock::getDescription() const {
std::stringstream s;
s << "Shock width: " << lShock / km << " km, ";
s << "Vup: " << vUp / km * sec << " km/s, ";
s << "Compression: " << compressionRatio;
return s.str();
}

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

OneDimensionalSphericalShock::OneDimensionalSphericalShock(double rShock, double vUp, double compressionRatio, double lShock, bool coolUpstream ){
setComp(compressionRatio);
setVup(vUp);
setShockwidth(lShock);
setShockRadius(rShock);
setCooling(coolUpstream);
}

Vector3d OneDimensionalSphericalShock::getField(const Vector3d &position) const {
double r = position.getR();
Vector3d e_r = position.getUnitVector();

double vDown = vUp / compressionRatio;
double a = (vUp + vDown) * 0.5;
double b = (vUp - vDown) * 0.5;

double v;
if (coolUpstream == true){

if (r <= rShock)
v = a - b * tanh((r-rShock) / lShock);
else
v = (a - b * tanh((r-rShock) / lShock)) * (rShock / r) * (rShock / r);

}
else
v = (a - b * tanh((r-rShock) / lShock)) * (rShock / r) * (rShock / r);

return v * e_r;
}

double OneDimensionalSphericalShock::getDivergence(const Vector3d &position) const {
double r = position.getR();

double vDown = vUp / compressionRatio;
double a = (vUp + vDown) * 0.5;
double b = (vUp - vDown) * 0.5;

double c = tanh((r-rShock) / lShock);

if (coolUpstream == true){
if (r <= rShock)
return 2 * a / r - 2 * b / r * c - b / lShock * (1 - c * c);
else
return -(rShock / r) * (rShock / r) * b / lShock * (1 - c * c);
}
else
return -(rShock / r) * (rShock / r) * b / lShock * (1 - c * c);

}

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

void OneDimensionalSphericalShock::setVup(double v) {
vUp = v;
return;
}
void OneDimensionalSphericalShock::setShockwidth(double w) {
lShock = w;
return;
}

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

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

double OneDimensionalSphericalShock::getComp() const {
return compressionRatio;
}

double OneDimensionalSphericalShock::getVup() const {
return vUp;
}

double OneDimensionalSphericalShock::getShockwidth() const {
return lShock;
}

double OneDimensionalSphericalShock::getShockRadius() const {
return rShock;
}

bool OneDimensionalSphericalShock::getCooling() const {
return coolUpstream;
}

std::string OneDimensionalSphericalShock::getDescription() const {
std::stringstream s;
s << "Shock width: " << lShock / km << " km, ";
s << "Shock radius: " << rShock / km << " km, ";
s << "Vup: " << vUp / km * sec << " km/s, ";
s << "Comp: " << compressionRatio;

return s.str();
}

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

ObliqueAdvectionShock::ObliqueAdvectionShock(double compressionRatio, double vXUp, double vY, double lShock) {
setComp(compressionRatio);
setVup(vXUp);
setVy(vY);
setShockwidth(lShock);
}

Vector3d ObliqueAdvectionShock::getField(const Vector3d &position) const {
double x = position.x;
double vXDown = vXUp / compressionRatio;

double a = (vXUp + vXDown) * 0.5;
double b = (vXUp - vXDown) * 0.5;

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

return v;
}

double ObliqueAdvectionShock::getDivergence(const Vector3d &position) const {
double x = position.x;
double vXDown = vXUp / compressionRatio;
// vy = const

double a = (vXUp + vXDown) * 0.5;
double b = (vXUp - vXDown) * 0.5;

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

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

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

void ObliqueAdvectionShock::setVy(double v) {
vY = v;
return;
}
void ObliqueAdvectionShock::setShockwidth(double w) {
lShock = w;
return;
}

double ObliqueAdvectionShock::getComp() const {
return compressionRatio;
}

double ObliqueAdvectionShock::getVup() const {
return vXUp;
}

double ObliqueAdvectionShock::getVy() const {
return vY;
}

double ObliqueAdvectionShock::getShockwidth() const {
return lShock;
}

std::string ObliqueAdvectionShock::getDescription() const {
std::stringstream s;
s << "Shock width: " << lShock / km << " km, ";
s << "Vx_up: " << vXUp / km * sec << " km/s, ";
s << "Vy: " << vY / km * sec << " km/s, ";
s << "Comp: " << compressionRatio;

return s.str();
}

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

Expand Down

0 comments on commit 026c6d4

Please sign in to comment.