Skip to content

Commit

Permalink
adlnTdlnP..
Browse files Browse the repository at this point in the history
  • Loading branch information
@jihenghu
jihenghu committed Mar 28, 2024
1 parent c99670e commit 46c1ba6
Show file tree
Hide file tree
Showing 4 changed files with 57 additions and 109 deletions.
16 changes: 8 additions & 8 deletions examples/2024-CLi-juno-mwr/run_juno_mwr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,20 +50,20 @@ def modify_atmosphere(
out.make_outputs(mesh, pin)


print(len(mesh.meshblocks()))
# print(len(mesh.meshblocks()))

mb = mesh.meshblock(0)


# for mb in mesh.meshblocks():
# print(mb.block_size.nx1)

# modify atmosphere
# dlnTdlnP=0
# pmin=0
# pmax=0
# adlnNH3dlnP=0
# mesh=modify_atmosphere(mesh,dlnTdlnP,pmin,pmax,"Temp")
# mesh=modify_atmosphere(mesh,adlnNH3dlnP,pmin,pmax,"NH3")
adlnNH3dlnP=30.0
pmin=1.E5
pmax=20.E5
mb.modify_dlnNH3dlnP(adlnNH3dlnP, pmin, pmax)

print(mb)


# run rt again
Expand Down
97 changes: 10 additions & 87 deletions python/pyathena.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,96 +78,23 @@ void init_athena(py::module &parent) {
}
});

// RegionSize
py::class_<RegionSize>(m, "RegionSize")
.def_property(
"x1min", [](RegionSize const &rs) { return rs.x1min; },
[](RegionSize &rs, Real x1min) { rs.x1min = x1min; })

.def_property(
"x2min", [](RegionSize const &rs) { return rs.x2min; },
[](RegionSize &rs, Real x2min) { rs.x2min = x2min; })

.def_property(
"x3min", [](RegionSize const &rs) { return rs.x3min; },
[](RegionSize &rs, Real x3min) { rs.x3min = x3min; })

.def_property(
"x1max", [](RegionSize const &rs) { return rs.x1max; },
[](RegionSize &rs, Real x1max) { rs.x1max = x1max; })

.def_property(
"x2max", [](RegionSize const &rs) { return rs.x2max; },
[](RegionSize &rs, Real x2max) { rs.x2max = x2max; })

.def_property(
"x3max", [](RegionSize const &rs) { return rs.x3max; },
[](RegionSize &rs, Real x3max) { rs.x3max = x3max; })

.def_property(
"nx1", [](RegionSize const &rs) { return rs.nx1; },
[](RegionSize &rs, int nx1) { rs.nx1 = nx1; })

.def_property(
"nx2", [](RegionSize const &rs) { return rs.nx2; },
[](RegionSize &rs, int nx2) { rs.nx2 = nx2; })

.def_property(
"nx3", [](RegionSize const &rs) { return rs.nx3; },
[](RegionSize &rs, int nx3) { rs.nx3 = nx3; })

.def_property(
"x1rat", [](RegionSize const &rs) { return rs.x1rat; },
[](RegionSize &rs, Real x1rat) { rs.x1rat = x1rat; })

.def_property(
"x2rat", [](RegionSize const &rs) { return rs.x2rat; },
[](RegionSize &rs, Real x2rat) { rs.x2rat = x2rat; })

.def_property(
"x3rat", [](RegionSize const &rs) { return rs.x3rat; },
[](RegionSize &rs, Real x3rat) { rs.x3rat = x3rat; });

// Mesh
py::class_<Mesh>(m, "Mesh")
.def(py::init<ParameterInput *, int>(), py::arg("pin"),
py::arg("mesh_only") = false)

.def("initialize",
[](Mesh &mesh, ParameterInput *pin) {
bool restart = false;

This comment has been minimized.

Sign in to view

Copy link
@jihenghu

jihenghu Mar 28, 2024

Author Collaborator

deleted by accident, recovered in next commit

// set up components
for (int b = 0; b < mesh.nblocal; ++b) {
MeshBlock *pmb = mesh.my_blocks(b);
pmb->pimpl = std::make_shared<MeshBlock::Impl>(pmb, pin);
}
mesh.Initialize(restart, pin);
})

.def("meshblocks",
[](Mesh &mesh) {
py::list lst;
for (size_t i = 0; i < mesh.nbtotal; ++i) {
lst.append(mesh.my_blocks(i));
}
return lst;
})
.def("initialize", [](Mesh &mesh, ParameterInput *pin) {
bool restart = false;

.def(
"meshblock",
[](Mesh &mesh, int n) {
if (n < 0 || n >= mesh.nbtotal) {
throw py::index_error();
}
return mesh.my_blocks(n);
},
py::return_value_policy::reference);
// set up components
for (int b = 0; b < mesh.nblocal; ++b) {
MeshBlock *pmb = mesh.my_blocks(b);
pmb->pimpl = std::make_shared<MeshBlock::Impl>(pmb, pin);
}
mesh.Initialize(restart, pin);
});

// MeshBlock
py::class_<MeshBlock>(m, "MeshBlock")
.def_readonly("block_size", &MeshBlock::block_size);

py::class_<MeshBlock>(m, "mesh_block")
.def("modify_dlnTdlnP", [](MeshBlock *mesh_block, Real adlnTdlnP, Real pmin, Real pmax) {
return modify_atmoshere_adlnTdlnP(mesh_block, adlnTdlnP, pmin, pmax);
Expand All @@ -177,11 +104,7 @@ void init_athena(py::module &parent) {
return modify_atmoshere_adlnNH3dlnP(mesh_block, adlnNH3dlnP, pmin, pmax);
});

//.def_readonly("inversion", [](MeshBlock const& pmb) {
// return pmb.pimpl->all_fits;
//});

// Outputs
// outputs
py::class_<Outputs>(m, "Outputs")
.def(py::init<Mesh *, ParameterInput *>(), py::arg("mesh"),
py::arg("pin"))
Expand Down
47 changes: 35 additions & 12 deletions src/utils/modify_atmoshere.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,34 +12,50 @@
#include <athena/hydro/hydro.hpp>
#include <athena/mesh/mesh.hpp>


// canoe
#include <air_parcel.hpp>
#include <constants.hpp>
#include <impl.hpp>
#include <dirty.hpp>

// snap
#include <snap/thermodynamics/thermodynamics.hpp>



// helper functions, will be moved in the future
int find_pressure_level_lesser_pybind(Real pres, AthenaArray<Real> const &w, int k,
int j, int is, int ie) {
for (int i = is; i <= ie; ++i)
if (w(IPR, k, j, i) < pres) return i;

return ie + 1;
}

// modify atmoshere with adlnTdlnP
void modify_atmoshere_adlnTdlnP(MeshBlock *pmb, Real adlnTdlnP, Real pmin, Real pmax) {
int is=pmb->is, js=pmb->js, ks=pmb->ks;
int ie=pmb->ie, je=pmb->je, ke=pmb->ke;

Hydro* phydro=pmb->phydro;
auto pthermo=Thermodynamics::GetInstance();
auto pcoord=pmb->pcoord;
Real H0 = pcoord->GetPressureScaleHeight();
Real dlnp = pcoord->dx1f(is) / H0;


// loop over all aircolumns
for (int k = ks; k <= ke; ++k) {
for (int j = js; j <= je; ++j) {
int ibegin = find_pressure_level_lesser(pmax, phydro->w, k, j, is, ie);
int iend = find_pressure_level_lesser(pmin, phydro->w, k, j, is, ie);
int ibegin = find_pressure_level_lesser_pybind(pmax, phydro->w, k, j, is, ie);
int iend = find_pressure_level_lesser_pybind(pmin, phydro->w, k, j, is, ie);

auto &&air = AirParcelHelper::gather_from_primitive(this, k, j, ibegin);
auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, ibegin);
air.ToMoleFraction();

for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat,
0., adlnTdlnP);
AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
AirParcelHelper::distribute_to_primitive(pmb, k, j, i + 1, air);
}
}
}
Expand All @@ -48,25 +64,32 @@ void modify_atmoshere_adlnTdlnP(MeshBlock *pmb, Real adlnTdlnP, Real pmin, Real
// modify atmoshere with adlnNH3dlnP
void modify_atmoshere_adlnNH3dlnP(MeshBlock *pmb, Real adlnNH3dlnP, Real pmin, Real pmax){

Hydro* phydro=pmb->phydro;
int is=pmb->is, js=pmb->js, ks=pmb->ks;
int ie=pmb->ie, je=pmb->je, ke=pmb->ke;

Hydro* phydro=pmb->phydro;
auto pthermo=Thermodynamics::GetInstance();
auto pcoord=pmb->pcoord;
Real H0 = pcoord->GetPressureScaleHeight();
Real dlnp = pcoord->dx1f(is) / H0;

// loop over all aircolumns
for (int k = ks; k <= ke; ++k) {
for (int j = js; j <= je; ++j) {
int ibegin = find_pressure_level_lesser(pmax, phydro->w, k, j, is, ie);
int iend = find_pressure_level_lesser(pmin, phydro->w, k, j, is, ie);
int ibegin = find_pressure_level_lesser_pybind(pmax, phydro->w, k, j, is, ie);
int iend = find_pressure_level_lesser_pybind(pmin, phydro->w, k, j, is, ie);

auto &&air = AirParcelHelper::gather_from_primitive(this, k, j, ibegin);
auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, ibegin);
air.ToMoleFraction();

for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat,
0., adlnNH3dlnP);
AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
AirParcelHelper::distribute_to_primitive(pmb, k, j, i + 1, air);
}
}
}
}
};



Expand Down
6 changes: 4 additions & 2 deletions src/utils/modify_atmoshere.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,14 @@
// athena
#include <athena/mesh/mesh.hpp>

// helper functions, will be moved in the future
int find_pressure_level_lesser_pybind(Real pres, AthenaArray<Real> const &w, int k, int j, int is, int ie) ;

// modify atmoshere with adlnTdlnP
void modify_atmoshere_adlnTdlnP(MeshBlock *pmb, Real adlnTdlnP, Real pmin, Real pmax);
void modify_atmoshere_adlnTdlnP(MeshBlock *pmb, Real adlnTdlnP, Real pmin, Real pmax) ;

// modify atmoshere with adlnNH3dlnP
void modify_atmoshere_adlnNH3dlnP(MeshBlock *pmb, Real adlnNH3dlnP, Real pmin, Real pmax);
void modify_atmoshere_adlnNH3dlnP(MeshBlock *pmb, Real adlnNH3dlnP, Real pmin, Real pmax) ;

#endif //SRC_UTILS_MODIFY_ATMOSPHERE_HPP_

0 comments on commit 46c1ba6

Please sign in to comment.