From 46c1ba68b411852fd8d736b64f9bdb0a079ecb1b Mon Sep 17 00:00:00 2001 From: jihenghu Date: Wed, 27 Mar 2024 21:53:23 -0400 Subject: [PATCH] adlnTdlnP.. --- examples/2024-CLi-juno-mwr/run_juno_mwr.py | 16 ++-- python/pyathena.cpp | 97 +++------------------- src/utils/modify_atmoshere.cpp | 47 ++++++++--- src/utils/modify_atmoshere.hpp | 6 +- 4 files changed, 57 insertions(+), 109 deletions(-) diff --git a/examples/2024-CLi-juno-mwr/run_juno_mwr.py b/examples/2024-CLi-juno-mwr/run_juno_mwr.py index f5aa9f42..263241e5 100644 --- a/examples/2024-CLi-juno-mwr/run_juno_mwr.py +++ b/examples/2024-CLi-juno-mwr/run_juno_mwr.py @@ -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 diff --git a/python/pyathena.cpp b/python/pyathena.cpp index 4c174910..fc99a6ca 100644 --- a/python/pyathena.cpp +++ b/python/pyathena.cpp @@ -78,96 +78,23 @@ void init_athena(py::module &parent) { } }); - // RegionSize - py::class_(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_(m, "Mesh") .def(py::init(), py::arg("pin"), py::arg("mesh_only") = false) - .def("initialize", - [](Mesh &mesh, ParameterInput *pin) { - bool restart = false; - - // set up components - for (int b = 0; b < mesh.nblocal; ++b) { - MeshBlock *pmb = mesh.my_blocks(b); - pmb->pimpl = std::make_shared(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(pmb, pin); + } + mesh.Initialize(restart, pin); + }); // MeshBlock - py::class_(m, "MeshBlock") - .def_readonly("block_size", &MeshBlock::block_size); - py::class_(m, "mesh_block") .def("modify_dlnTdlnP", [](MeshBlock *mesh_block, Real adlnTdlnP, Real pmin, Real pmax) { return modify_atmoshere_adlnTdlnP(mesh_block, adlnTdlnP, pmin, pmax); @@ -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_(m, "Outputs") .def(py::init(), py::arg("mesh"), py::arg("pin")) diff --git a/src/utils/modify_atmoshere.cpp b/src/utils/modify_atmoshere.cpp index 6695913a..57f80819 100644 --- a/src/utils/modify_atmoshere.cpp +++ b/src/utils/modify_atmoshere.cpp @@ -12,34 +12,50 @@ #include #include - // canoe #include #include #include -#include // snap #include + + +// helper functions, will be moved in the future +int find_pressure_level_lesser_pybind(Real pres, AthenaArray 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); } } } @@ -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); } } } -} +}; diff --git a/src/utils/modify_atmoshere.hpp b/src/utils/modify_atmoshere.hpp index 860bdf90..b71d1940 100644 --- a/src/utils/modify_atmoshere.hpp +++ b/src/utils/modify_atmoshere.hpp @@ -11,12 +11,14 @@ // athena #include +// helper functions, will be moved in the future +int find_pressure_level_lesser_pybind(Real pres, AthenaArray 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_