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.

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

lavaplanet input and cpp files #115

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
Changes from 1 commit
Commits
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
Prev Previous commit
Next Next commit
add SiO svp function
  • Loading branch information
gitlinffff committed Jan 22, 2024
commit 64fc11b9c8060d460ce519b439c1fc3acacb8084
366 changes: 202 additions & 164 deletions examples/2023-Linf-lava/lava.cpp
Original file line number Diff line number Diff line change
@@ -33,180 +33,213 @@
// utils
#include <utils/fileio.hpp>

//Real grav, P0, T0, Tmin;
//int iH2O;
//
//
//
//
//void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
// AllocateUserOutputVariables(4 + NVAPOR);
// SetUserOutputVariableName(0, "temp");
// SetUserOutputVariableName(1, "theta");
// SetUserOutputVariableName(2, "thetav");
// SetUserOutputVariableName(3, "mse");
// for (int n = 1; n <= NVAPOR; ++n) {
// std::string name = "rh" + std::to_string(n);
// SetUserOutputVariableName(3 + n, name.c_str());
// }
//}
//
//void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
// auto pthermo = Thermodynamics::GetInstance();
//
// for (int k = ks; k <= ke; ++k)
// for (int j = js; j <= je; ++j)
// for (int i = is; i <= ie; ++i) {
// user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i);
// user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i);
// // theta_v
// user_out_var(2, k, j, i) =
// user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i);
// // mse
// user_out_var(3, k, j, i) =
// pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
// for (int n = 1; n <= NVAPOR; ++n)
// user_out_var(3 + n, k, j, i) =
// pthermo->RelativeHumidity(this, n, k, j, i);
// }
//}
//
//void Forcing(MeshBlock *pmb, Real const time, Real const dt,
// AthenaArray<Real> const &w, AthenaArray<Real> const &r,
// AthenaArray<Real> const &bcc, AthenaArray<Real> &du,
// AthenaArray<Real> &s) {
// int is = pmb->is, js = pmb->js, ks = pmb->ks;
// int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
// auto pthermo = Thermodynamics::GetInstance();
//
// for (int k = ks; k <= ke; ++k)
// for (int j = js; j <= je; ++j)
// for (int i = is; i <= ie; ++i) {
// auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
//
// Real cv = pthermo->GetCvMass(air, 0);
//
// // if (w(IPR, k, j, i) < prad) {
// // du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv *
// // (1. + 1.E-4 * sin(2. * M_PI * rand() /
// // RAND_MAX));
// // }
//
// // if (air.w[IDN] < Tmin) {
// // u(IEN,k,j,i) += w(IDN,k,j,i)*cv*(Tmin - temp)/sponge_tau*dt;
// // }
// }
//}
//
//void Mesh::InitUserMeshData(ParameterInput *pin) {
// grav = -pin->GetReal("hydro", "grav_acc1");
//
// P0 = pin->GetReal("problem", "P0");
// T0 = pin->GetReal("problem", "T0");
//
// Tmin = pin->GetReal("problem", "Tmin");
//
// // index
// auto pindex = IndexMap::GetInstance();
// iH2O = pindex->GetVaporId("H2O");
// EnrollUserExplicitSourceFunction(Forcing);
//}
Real grav, P0, T0, Tmin;
int iSiO = 2; // assign an VaporID to iSiO (should have got from pindex)
int iH2O;




void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
AllocateUserOutputVariables(4 + NVAPOR);
SetUserOutputVariableName(0, "temp");
SetUserOutputVariableName(1, "theta");
SetUserOutputVariableName(2, "thetav");
SetUserOutputVariableName(3, "mse");
for (int n = 1; n <= NVAPOR; ++n) {
std::string name = "rh" + std::to_string(n);
SetUserOutputVariableName(3 + n, name.c_str());
}
}

void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
auto pthermo = Thermodynamics::GetInstance();

for (int k = ks; k <= ke; ++k)
for (int j = js; j <= je; ++j)
for (int i = is; i <= ie; ++i) {
user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i);
user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i);
// theta_v
user_out_var(2, k, j, i) =
user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i);
// mse
user_out_var(3, k, j, i) =
pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
for (int n = 1; n <= NVAPOR; ++n)
user_out_var(3 + n, k, j, i) =
pthermo->RelativeHumidity(this, n, k, j, i);
}
}

void Forcing(MeshBlock *pmb, Real const time, Real const dt,
AthenaArray<Real> const &w, AthenaArray<Real> const &r,
AthenaArray<Real> const &bcc, AthenaArray<Real> &du,
AthenaArray<Real> &s) {
int is = pmb->is, js = pmb->js, ks = pmb->ks;
int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
auto pthermo = Thermodynamics::GetInstance();

for (int k = ks; k <= ke; ++k)
for (int j = js; j <= je; ++j)
for (int i = is; i <= ie; ++i) {
auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);

Real cv = pthermo->GetCvMass(air, 0);

// if (w(IPR, k, j, i) < prad) {
// du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv *
// (1. + 1.E-4 * sin(2. * M_PI * rand() /
// RAND_MAX));
// }

// if (air.w[IDN] < Tmin) {
// u(IEN,k,j,i) += w(IDN,k,j,i)*cv*(Tmin - temp)/sponge_tau*dt;
// }
}
}

void Mesh::InitUserMeshData(ParameterInput *pin) {
grav = -pin->GetReal("hydro", "grav_acc1");

P0 = pin->GetReal("problem", "P0");
T0 = pin->GetReal("problem", "T0");

Tmin = pin->GetReal("problem", "Tmin");

// index
auto pindex = IndexMap::GetInstance();
iH2O = pindex->GetVaporId("H2O");
EnrollUserExplicitSourceFunction(Forcing);
}



void MeshBlock::ProblemGenerator(ParameterInput *pin) {
srand(Globals::my_rank + time(0));

Application::Logger app("main");
app->Log("ProblemGenerator: jupiter_crm");

auto pthermo = Thermodynamics::GetInstance();

// mesh limits
Real x1min = pmy_mesh->mesh_size.x1min;
Real x1max = pmy_mesh->mesh_size.x1max;

// request temperature and pressure
app->Log("request T", T0);
app->Log("request P", P0);

void MeshBlock::ProblemGenerator(ParameterInput *pin){ // sets initial conditions
// thermodynamic constants
Real gamma = pin->GetReal("hydro", "gamma");
Real Rd = pthermo->GetRd();
Real cp = gamma / (gamma - 1.) * Rd;


std::string input_atm_path = "/home/linfel/canoe_1/examples/2023-Linf-lava/balance_atm.txt";
DataVector atm_data = read_data_vector(input_atm_path);



///////////////////////////////////////////////////////////////////////////////////////////////
// set up an adiabatic atmosphere
//int max_iter = 400, iter = 0;
//Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01);

// Access the data using keys
for (const auto& entry : atm_data) {
const std::string& key = entry.first;
const std::vector<double>& values = entry.second;
AirParcel air(AirParcel::Type::MoleFrac);

// Print the key (column header)
std::cout << key << ": ";
// estimate surface temperature and pressure
//Real Ts = T0 - grav / cp * x1min;
//Real Ps = P0 * pow(Ts / T0, cp / Rd);
//Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6;

// Print the values in the vector
for (double value : values) {
std::cout << value << " ";
}
//while (iter++ < max_iter) {
// // read in vapors
// air.w[iH2O] = xH2O;
// air.w[IPR] = Ps;
// air.w[IDN] = Ts;

// // stop at just above P0
// for (int i = is; i <= ie; ++i) {
// pthermo->Extrapolate(&air, pcoord->dx1f(i),
// Thermodynamics::Method::PseudoAdiabat, grav);
// if (air.w[IPR] < P0) break;
// }

// // make up for the difference
// Ts += T0 - air.w[IDN];
// if (std::abs(T0 - air.w[IDN]) < Ttol) break;

// app->Log("Iteration #", iter);
// app->Log("T", air.w[IDN]);
//}

//if (iter > max_iter) {
// throw RuntimeError("ProblemGenerator", "maximum iteration reached");
//}
//////////////////////////////////////////////////////////////////////////////////////


std::cout << std::endl;
Real dz, **w1, *z1;
int nx1 = 80; // int or size_t ? The function 'interpn' require a size_t type
int nfield = 3; // number of fields (PRE,TEM,SiO)
z1 = new Real [nx1];
//NewCArray(w1, nx1, nfield); // which header should be included for 'NewCArray'?
w1 = new Real *[nx1]; // Allocate memory for w1
for (int i = 0; i < nx1; ++i) {
w1[i] = new Real[nfield];
}
}
//std::cout << " IPR " << IPR << " IDN " << IDN << " iSiO " << iSiO << std::endl;
for (int i = 0; i < nx1; ++i){
z1[i] = atm_data["HGT"][i];
w1[i][IPR] = atm_data["PRE"][i];
w1[i][IDN] = atm_data["TEM"][i];
w1[i][iSiO] = atm_data["SiO"][i];
//std::cout << "i=" << i << " HGT " << z1[i] << " PRE " << w1[i][IPR] << " TEM " << w1[i][IDN] << " SIO " << w1[i][iSiO] << std::endl;
}

// construct atmosphere from bottom up
for (int k = ks; k <= ke; ++k){
for (int j = js; j <= je; ++j){
// half a grid to cell center
pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
Thermodynamics::Method::ReversibleAdiabat, grav);


for (int i = is; i <= ie; ++i) {
// interpolation to coord grid
Real buf[nfield];
interpn(buf, &pcoord->x1v(i), *w1, z1, &nx1, 1, nfield);
//std::cout << "i=" << i << ", j=" << j << ", k=" << k << " buf0 " << buf[0] << " buf1 " << buf[1] << " buf2 " << buf[2] << std::endl;
// set physical state for an air parcel
// air.SetZero();
// air.w[IPR] = buf[];
// air.w[IDN] = buf[];
// air.w[iSiO] = buf[];
//
// //std::cout << "i=" << i << ", j=" << j << ", k=" << k << std::endl;
// //air.w[IVX] = 1. * sin(2. * M_PI * rand() / RAND_MAX);
// AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
// pthermo->Extrapolate(&air, pcoord->dx1f(i),
// Thermodynamics::Method::PseudoAdiabat, grav,
// 1.e-5);
}
}
}











//void MeshBlock::ProblemGenerator(ParameterInput *pin) {
// srand(Globals::my_rank + time(0));
//
// Application::Logger app("main");
// app->Log("ProblemGenerator: jupiter_crm");
//
// auto pthermo = Thermodynamics::GetInstance();
//
// // mesh limits
// Real x1min = pmy_mesh->mesh_size.x1min;
// Real x1max = pmy_mesh->mesh_size.x1max;
//
// // request temperature and pressure
// app->Log("request T", T0);
// app->Log("request P", P0);
//
// // thermodynamic constants
// Real gamma = pin->GetReal("hydro", "gamma");
// Real Rd = pthermo->GetRd();
// Real cp = gamma / (gamma - 1.) * Rd;
//
// // set up an adiabatic atmosphere
// int max_iter = 400, iter = 0;
// Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01);
//
// AirParcel air(AirParcel::Type::MoleFrac);
//
// // estimate surface temperature and pressure
// Real Ts = T0 - grav / cp * x1min;
// Real Ps = P0 * pow(Ts / T0, cp / Rd);
// Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6;
//
// while (iter++ < max_iter) {
// // read in vapors
// air.w[iH2O] = xH2O;
// air.w[IPR] = Ps;
// air.w[IDN] = Ts;
//
// // stop at just above P0
// for (int i = is; i <= ie; ++i) {
// pthermo->Extrapolate(&air, pcoord->dx1f(i),
// Thermodynamics::Method::PseudoAdiabat, grav);
// if (air.w[IPR] < P0) break;
// }
//
// // make up for the difference
// Ts += T0 - air.w[IDN];
// if (std::abs(T0 - air.w[IDN]) < Ttol) break;
//
// app->Log("Iteration #", iter);
// app->Log("T", air.w[IDN]);
// }
//
// if (iter > max_iter) {
// throw RuntimeError("ProblemGenerator", "maximum iteration reached");
// }
//
// // construct atmosphere from bottom up
// for (int k = ks; k <= ke; ++k)
// for (int j = js; j <= je; ++j) {
// air.SetZero();
// air.w[iH2O] = xH2O;
// air.w[IPR] = Ps;
// air.w[IDN] = Ts;
//
// // half a grid to cell center
// pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
// Thermodynamics::Method::ReversibleAdiabat, grav);
//
// int i = is;
// for (; i <= ie; ++i) {
// if (air.w[IDN] < Tmin) break;
// air.w[IVX] = 1. * sin(2. * M_PI * rand() / RAND_MAX);
// air.w[IVX] = 1. * sin(2. * M_PI * rand() / RAND_MAX);
// AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
// pthermo->Extrapolate(&air, pcoord->dx1f(i),
// Thermodynamics::Method::PseudoAdiabat, grav,
@@ -226,4 +259,9 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin){ // sets initial conditio
//
// pimpl->prad->CalFlux(this, k, j, is, ie + 1);
// }
//}
delete[] z1;
for (int i = 0; i < nx1; ++i) {
delete[] w1[i];
}
delete[] w1;
}
Loading