Skip to content

Commit

Permalink
Merge pull request fmihpc#960 from vetarvus/KHB
Browse files Browse the repository at this point in the history 
Updates to the Kelvin-Helmholtz instability project
  • Loading branch information
@markusbattarbee
markusbattarbee authored May 13, 2024
2 parents e2d9b46 + 2fea9d2 commit da5681a
Show file tree
Hide file tree
Showing 8 changed files with 257 additions and 48 deletions.
34 changes: 17 additions & 17 deletions projects/KHB/KHB.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,14 +15,14 @@ system_write_distribution_zline_stride = 0

[gridbuilder]
x_length = 20
y_length = 1
z_length = 20
y_length = 20
z_length = 1
x_min = -5.0e6
x_max = 5.0e6
y_min = 0.0
y_max = 5.0e5
z_min = 0.0
z_max = 1.0e7
y_min = 0.0
y_max = 1.0e7
z_min = 0.0
z_max = 5.0e5
vx_min = -600000.0
vx_max = +600000.0
vy_min = -600000.0
Expand Down Expand Up @@ -70,25 +70,25 @@ diagnostic = populations_vg_blocks
minValue = 1.0e-16

[KHB]
P = 2.765276873575132E-10

Vx1 = 0.0
Vy1 = 0.0
Vz1 = 4.0e5
Vy1 = 4.0e5
Vz1 = 0.0
Bx1 = 0.0
By1 = 1.0e-9
Bz1 = 0.0
T1 = 1.0e5
By1 = 0.0
Bz1 = 1.0e-9
rho1 = 1.0e6

Vx2 = 0.0
Vy2 = 0.0
Vz2 = -4.0e5
Vy2 = -4.0e5
Vz2 = 0.0
Bx2 = 0.0
By2 = 1.0e-9
Bz2 = 0.0
T2 = 1.0e5
By2 = 0.0
Bz2 = 1.0e-9
rho2 = 2.0e6

lambda = 3.3333333e6
lambda = 4.0e3
amp = 0.0
offset = 0.0
transitionWidth = 1.0e6
Expand Down
92 changes: 68 additions & 24 deletions projects/KHB/KHB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,17 +33,16 @@

namespace projects {
using namespace std;
KHB::KHB(): Project() { }
KHB::KHB(): TriAxisSearch() { }
KHB::~KHB() { }

bool KHB::initialize(void) {return Project::initialize();}

void KHB::addParameters() {
typedef Readparameters RP;
RP::add("KHB.P", "Constant total pressure (thermal+magnetic), used to determine the temperature profile (Pa)", 0.0);
RP::add("KHB.rho1", "Number density, this->TOP state (m^-3)", 0.0);
RP::add("KHB.rho2", "Number density, this->BOTTOM state (m^-3)", 0.0);
RP::add("KHB.T1", "Temperature, this->TOP state (K)", 0.0);
RP::add("KHB.T2", "Temperature, this->BOTTOM state (K)", 0.0);
RP::add("KHB.Vx1", "Bulk velocity x component, this->TOP state (m/s)", 0.0);
RP::add("KHB.Vx2", "Bulk velocity x component, this->BOTTOM state (m/s)", 0.0);
RP::add("KHB.Vy1", "Bulk velocity y component, this->TOP state (m/s)", 0.0);
Expand All @@ -57,9 +56,11 @@ namespace projects {
RP::add("KHB.Bz1", "Magnetic field z component, this->TOP state (T)", 0.0);
RP::add("KHB.Bz2", "Magnetic field z component, this->BOTTOM state (T)", 0.0);
RP::add("KHB.lambda", "Initial perturbation wavelength (m)", 0.0);
RP::add("KHB.amp", "Initial perturbation amplitude (m)", 0.0);
RP::add("KHB.amp", "Initial velocity perturbation amplitude (m s^-1)", 0.0);
RP::add("KHB.offset", "Boundaries offset from 0 (m)", 0.0);
RP::add("KHB.transitionWidth", "Width of tanh transition for all changing values", 0.0);
RP::add("KHB.harmonics", "Number of harmonics of lambda included in the initial perturbation", 0);
RP::add("KHB.randomPhase", "If true, set a random phase for each mode of the initial perturbation. Seed set via project_common.seed", 0);
}

void KHB::getParameters() {
Expand All @@ -71,10 +72,9 @@ namespace projects {
abort();
}

RP::get("KHB.P", this->P);
RP::get("KHB.rho1", this->rho[this->TOP]);
RP::get("KHB.rho2", this->rho[this->BOTTOM]);
RP::get("KHB.T1", this->T[this->TOP]);
RP::get("KHB.T2", this->T[this->BOTTOM]);
RP::get("KHB.Vx1", this->Vx[this->TOP]);
RP::get("KHB.Vx2", this->Vx[this->BOTTOM]);
RP::get("KHB.Vy1", this->Vy[this->TOP]);
Expand All @@ -91,37 +91,81 @@ namespace projects {
RP::get("KHB.amp", this->amp);
RP::get("KHB.offset", this->offset);
RP::get("KHB.transitionWidth", this->transitionWidth);
RP::get("KHB.harmonics", this->harmonics);
RP::get("KHB.randomPhase", this->randomPhase);
}


Real KHB::profile(creal top, creal bottom, creal x, creal z) const {
Real KHB::profile(creal top, creal bottom, creal x) const {
if(top == bottom) {
return top;
}
if(this->offset != 0.0) {
return 0.5 * ((top-bottom) * (
tanh((x + this->offset + this->amp * cos(2.0*M_PI*z/this->lambda))/this->transitionWidth) -
tanh((x-(this->offset + this->amp * cos(2.0*M_PI*z/this->lambda)))/this->transitionWidth) -1) + top+bottom);
tanh((x + this->offset)/this->transitionWidth) -
tanh((x - this->offset)/this->transitionWidth) -1) + top+bottom);
} else {
return 0.5 * ((top-bottom) * tanh(x/this->transitionWidth) + top+bottom);
}
}

Real KHB::getDistribValue(creal& x, creal& z, creal& vx, creal& vy, creal& vz, const uint popID) const {
creal mass = physicalconstants::MASS_PROTON;
creal kb = physicalconstants::K_B;
Real rho = profile(this->rho[this->BOTTOM], this->rho[this->TOP], x, z);
Real T = profile(this->T[this->BOTTOM], this->T[this->TOP], x, z);
Real Vx = profile(this->Vx[this->BOTTOM], this->Vx[this->TOP], x, z);
Real Vy = profile(this->Vy[this->BOTTOM], this->Vy[this->TOP], x, z);
Real Vz = profile(this->Vz[this->BOTTOM], this->Vz[this->TOP], x, z);

return rho * pow(mass / (2.0 * M_PI * kb * T), 1.5) *
exp(- mass * (pow(vx - Vx, 2.0) + pow(vy - Vy, 2.0) + pow(vz - Vz, 2.0)) / (2.0 * kb * T));
inline vector<std::array<Real, 3> > KHB::getV0(
creal x,
creal y,
creal z,
const uint popID
) const {
Real Vx = profile(this->Vx[this->BOTTOM], this->Vx[this->TOP], x);
Real Vy = profile(this->Vy[this->BOTTOM], this->Vy[this->TOP], x);
Real Vz = profile(this->Vz[this->BOTTOM], this->Vz[this->TOP], x);
(void)y,z,popID;

// add an initial velocity perturbation to Vx
// initialize RNG for calculating random phases for the initial perturbation
std::default_random_engine rndState;
setRandomSeed(0,rndState);
Real phase = 0.0;

// add each mode to the initial perturbation
for (int i=0; i<=this->harmonics; i++) {
if (this->randomPhase) {
phase = 2.0 * M_PI * getRandomNumber(rndState);
}

if (this->offset != 0.0) {
Vx += this->amp * sin(2.0 * (i + 1) * M_PI * y / this->lambda + phase) * (exp(-pow((x + this->offset) / this->transitionWidth,2)) + exp(-pow((x - this->offset) / this->transitionWidth,2)));
} else {
Vx += this->amp * sin(2.0 * (i + 1) * M_PI * y / this->lambda + phase) * exp(-pow(x / this->transitionWidth,2));
}
}

vector<std::array<Real, 3> > centerPoints;
std::array<Real, 3> V0 {{Vx,Vy,Vz}};
centerPoints.push_back(V0);
return centerPoints;
}

inline Real KHB::getDistribValue(creal& x, creal& y, creal& z, creal& vx, creal& vy, creal& vz, const uint popID) const {
creal mass = getObjectWrapper().particleSpecies[popID].mass;
Real rho = profile(this->rho[this->BOTTOM], this->rho[this->TOP], x);
std::array<Real, 3> initV0 = this->getV0(x, y, z, popID)[0];
Real Vx = initV0[0];
Real Vy = initV0[1];
Real Vz = initV0[2];

// calculate the temperature such that the total pressure is constant across the domain
Real mu0 = physicalconstants::MU_0;
Real Bx = profile(this->Bx[this->BOTTOM], this->Bx[this->TOP], x);
Real By = profile(this->By[this->BOTTOM], this->By[this->TOP], x);
Real Bz = profile(this->Bz[this->BOTTOM], this->Bz[this->TOP], x);
Real kbT = (this->P - 0.5 * (Bx * Bx + By * By + Bz * Bz) / mu0) / rho;

return rho * pow(mass / (2.0 * M_PI * kbT), 1.5) *
exp(- mass * (pow(vx - Vx, 2.0) + pow(vy - Vy, 2.0) + pow(vz - Vz, 2.0)) / (2.0 * kbT));
}

Real KHB::calcPhaseSpaceDensity(creal& x, creal& y, creal& z, creal& dx, creal& dy, creal& dz, creal& vx, creal& vy, creal& vz, creal& dvx, creal& dvy, creal& dvz,const uint popID) const {
return getDistribValue(x+0.5*dx, z+0.5*dz, vx+0.5*dvx, vy+0.5*dvy, vz+0.5*dvz, popID);
return getDistribValue(x+0.5*dx, y+0.5*dy, z+0.5*dz, vx+0.5*dvx, vy+0.5*dvy, vz+0.5*dvz, popID);
}

void KHB::calcCellParameters(spatial_cell::SpatialCell* cell,creal& t) { }
Expand All @@ -143,9 +187,9 @@ namespace projects {
const std::array<Real, 3> xyz = perBGrid.getPhysicalCoords(x, y, z);
std::array<Real, fsgrids::bfield::N_BFIELD>* cell = perBGrid.get(x, y, z);

cell->at(fsgrids::bfield::PERBX) = profile(this->Bx[this->BOTTOM], this->Bx[this->TOP], xyz[0]+0.5*perBGrid.DX, xyz[2]+0.5*perBGrid.DZ);
cell->at(fsgrids::bfield::PERBY) = profile(this->By[this->BOTTOM], this->By[this->TOP], xyz[0]+0.5*perBGrid.DX, xyz[2]+0.5*perBGrid.DZ);
cell->at(fsgrids::bfield::PERBZ) = profile(this->Bz[this->BOTTOM], this->Bz[this->TOP], xyz[0]+0.5*perBGrid.DX, xyz[2]+0.5*perBGrid.DZ);
cell->at(fsgrids::bfield::PERBX) = profile(this->Bx[this->BOTTOM], this->Bx[this->TOP], xyz[0]+0.5*perBGrid.DX);
cell->at(fsgrids::bfield::PERBY) = profile(this->By[this->BOTTOM], this->By[this->TOP], xyz[0]+0.5*perBGrid.DX);
cell->at(fsgrids::bfield::PERBZ) = profile(this->Bz[this->BOTTOM], this->Bz[this->TOP], xyz[0]+0.5*perBGrid.DX);
}
}
}
Expand Down
18 changes: 13 additions & 5 deletions projects/KHB/KHB.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,10 @@
#include <stdlib.h>

#include "../../definitions.h"
#include "../project.h"
#include "../projectTriAxisSearch.h"

namespace projects {
class KHB: public Project {
class KHB: public TriAxisSearch {
public:
KHB();
virtual ~KHB();
Expand All @@ -38,6 +38,12 @@ namespace projects {
static void addParameters(void);
virtual void getParameters(void);
virtual void calcCellParameters(spatial_cell::SpatialCell* cell,creal& t);
virtual std::vector<std::array<Real, 3> > getV0(
creal x,
creal y,
creal z,
const uint popID
) const;
virtual Real calcPhaseSpaceDensity(
creal& x, creal& y, creal& z,
creal& dx, creal& dy, creal& dz,
Expand All @@ -52,17 +58,17 @@ namespace projects {
);
protected:
Real getDistribValue(
creal& x, creal& z,
creal& x, creal& y, creal& z,
creal& vx, creal& vy, creal& vz,
const uint popID) const;
Real profile(creal top, creal bottom, creal x, creal z) const;
Real profile(creal top, creal bottom, creal x) const;

enum {
TOP,
BOTTOM
};
Real P;
Real rho[2];
Real T[2];
Real Vx[2];
Real Vy[2];
Real Vz[2];
Expand All @@ -73,6 +79,8 @@ namespace projects {
Real amp;
Real offset;
Real transitionWidth;
int harmonics;
bool randomPhase;
}; // class KHB
} // namespace
#endif
2 changes: 1 addition & 1 deletion projects/KHB/mxm.dat
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
0.0 1.0e6 1.0e7 0.0 0.0 4.0e6 1.0e-11 1.0e-7 1.0e-11
0.0 1.0e6 2.0e7 0.0 4.0e6 0.0 0.0 0.0 1.0e-9
2 changes: 1 addition & 1 deletion projects/KHB/mxp.dat
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
0.0 4.0e6 1.0e7 0.0 0.0 -4.0e6 1.0e-11 1.0e-7 1.0e-11
0.0 2.0e6 1.0e7 0.0 -4.0e6 0.0 0.0 0.0 1.0e-9
Loading

0 comments on commit da5681a

Please sign in to comment.