本节解读一下juno_mwr_eq示例的具体实现和实验使用到的理论和假设。
情景描述
Juno Mission 朱诺（Juno）是一项由美国宇航局（NASA）进行的木星探测任务, 于2011年8月5日从佛罗里达州的卡纳维拉尔角空军基地发射,于2016年7月4日成功进入了木星轨道, 主要目标是深入研究木星的内部结构、大气层和磁场，以帮助科学家们了解这个巨大气体行星的形成和演化过程。 朱诺搭载了多种科学仪器，包括微波辐射计(MWR)、微波雷达、飞行器磁力计、高能粒子探测器等，用于测量木星大气的温度、气体成分、磁场和辐射环境等参数。 极化轨道：
朱诺以一个相对较小的轨道倾角环绕木星，被称为“极化轨道”，以减少朱诺受到木星辐射带的辐射影响，并允许更全面地研究木星的极区域。
Juno的轨道离心率较大，可以使其再近木点钻进木星辐射带进行大气辐射的观测，减少木星磁场的干扰。
已知木星大气特性：伽利略号探测器对木星的穿透性探测得到了一系列的probe in-situ测了结果，其中包括，大气气体成分主要是水和氨气；1 bar 高度的气温约为169 K。 模拟需要实现的步骤：
- 输入已知的探测结果和热力学参数；
- 基于绝热假设构建大气温湿度廓线和气体含量廓线；
- 辐射传输模拟，计算天顶微波亮温。
运行个例
构建项目和编译
$ cd
(env) jihenghu@canopy:~ [11:31]$ mkdir -p Juno_MWR_Proj/build
(env) jihenghu@canopy:~ [11:31]$ cd Juno_MWR_Proj/build/
(env) jihenghu@canopy:build [11:32]$ cmake ~/canoe/ -DTASK=juno .
(env) jihenghu@canopy:build [11:32]$ make -j8 ## 8核编译
准备所需要的文件：
$ cd ..
(env) jihenghu@canopy:Juno_MWR_Proj [11:35]$ ln -s build/bin/juno_mwr.release .
(env) jihenghu@canopy:Juno_MWR_Proj [11:35]$ cp build/bin/juno_mwr.inp .
(env) jihenghu@canopy:Juno_MWR_Proj [11:36]$ cp build/bin/mwr_channels.yaml .
输入和配置文件juno_mwr.inp | mwr_channels.yaml我们会在后面介绍其具体意义。
运行
(env) jihenghu@canopy:Juno_MWR_Proj [11:40]$ ./juno_mwr.release -i juno_mwr.inp
######################################################
## CANOE MODEL STARTS ##
######################################################
Log, "2023-11-15 11:40:57", canoe, 2., "Installing monitor canoe"
Log, "2023-11-15 11:40:57", canoe, 2.1., "Initialize IndexMap"
......
Log, "2023-11-15 11:40:59", snap, 47.1., "Destroy Thermodynamics"
Log, "2023-11-15 11:40:59", canoe, 48.1., "Destroy IndexMap"
Log, "2023-11-15 11:40:59", outputs, 49.1., "Destroy MetadataTable"
Terminating on success
cpu time used = 1.6321420000000000e+00 (s)
检查输出：
$ tree
.
├── build/
├── juno_mwr.inp
├── juno_mwr.out2.00000.nc
├── juno_mwr.out3.00000.nc
├── juno_mwr.out4.00000.nc
├── juno_mwr.out5.00000.nc
├── juno_mwr.release -> build/bin/juno_mwr.release
└── mwr_channels.yaml
$ ln -s build/bin/combine.py
$ python3 combine.py -o test. ## 将输出合并成一个nc文件
$ tree
├── build/
├── combine.py -> build/bin/combine.py
├── combine_rules
├── juno_mwr.inp
├── juno_mwr.release -> build/bin/juno_mwr.release
├── juno_mwr-test-main.nc
└── mwr_channels.yam
结果分析
- 大气廓线输出：包括温度、位温、两种气体的浓度廓线。
- 辐射传输模拟：包括6个通道的天顶亮温，临边昏暗和权重廓线
Note:这里的权重指的是：
$w(z)=-e^{- \tau } d\tau /dz $
源码分析
比较特殊的头文件如下：
#include <snap/thermodynamics/vapors/sodium_vapors.hpp>
#include <astro/Jupiter/jup_fletcher16_cirs.hpp>
#include <astro/planets.hpp>
//radiative transfer
#include <harp/radiation.hpp>
#include <harp/radiation_band.hpp>
// tracer
#include <tracer/tracer.hpp>
// opacity
#include <opacity/Giants/microwave/mwr_absorbers.hpp>
// inversion
#include <inversion/profile_inversion.hpp>
// special includes
#include "juno_mwr_specs.hpp"
全局变量
Real grav, P0, T0, Tmin, clat;
Real xHe, xCH4, xH2S, xNa, xKCl, metallicity; // 痕量气体的含量
设定输出的诊断量，常规的温度和能量，加上几种condensible气体的相对湿度。这里只有H2O和NH3两种（NHYDRO=2）。
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());
}
}
NVAPOR,NHYDRO的定义在defs.hpp.in
输出前的计算
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 Mesh::InitUserMeshData(ParameterInput *pin) {
grav = -pin->GetReal("hydro", "grav_acc1");
P0 = pin->GetReal("mesh", "ReferencePressure");
T0 = pin->GetReal("problem", "T1bar");
Tmin = pin->GetReal("problem", "Tmin");
clat = pin->GetOrAddReal("problem", "clat", 0.);
xH2S = pin->GetReal("problem", "xH2S");
metallicity = pin->GetOrAddReal("problem", "metallicity", 0.);
xNa = pin->GetReal("problem", "xNa");
xNa *= pow(10., metallicity);
xKCl = pin->GetReal("problem", "xKCl");
xKCl *= pow(10., metallicity);
xHe = pin->GetReal("problem", "xHe");
xCH4 = pin->GetReal("problem", "xCH4");
}
ProblemGenerator
这个个例比较复杂，我们逐块解释：
- 数据读取
Application::Logger app("main"); // 日志注册
app->Log("ProblemGenerator: juno");
auto pthermo = Thermodynamics::GetInstance(); // 获取热力学指针对象
// mesh limits
Real x1min = pmy_mesh->mesh_size.x1min;
Real x1max = pmy_mesh->mesh_size.x1max;
Real H0 = pcoord->GetPressureScaleHeight(); // 大气标高H0
这里读取了模拟的垂直坐标的范围和大气标高H0： Input文件的范围是 x1min = -270.E3 到 x1max = 200.E3 的高度范围。 这里说明一下模拟使用的坐标系统： \(\begin{align} z^*= H_0 \cdot ln\frac{P}{P_0} ; P_0= 1 bar \end{align}\)
- 对于等温大气：$z=z^*=H_0 \cdot ln\frac{P}{P_0}$；
- 对于等熵大气：$z !=z^*$
// 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;
// index
auto pindex = IndexMap::GetInstance();
int iH2O = pindex->GetVaporId("H2O");
int iNH3 = pindex->GetVaporId("NH3");
app->Log("index of H2O", iH2O);
app->Log("index of NH3", iNH3);
迭代求算$T_s$, $P_s$:
- ①计算地表气压和地表气温（干绝热模型）： $ P_s=P_0 \cdot e^{-x1min/H0} $；$ T_s=T_0 \cdot \frac{P_s}{P_0}^{R_d/C_p} $ ;
- ②从地表开始构建一条绝热廓线到1bar的位置（湿绝热），获得是$T_0’$ ;
- 计算1bar处的温度偏差$\Delta T=T_0’ - T_0$，修正Ts：$ T_s’=T_s- \Delta T $，
- ③自下而上再次构建湿绝热廓线，获得新的T_0”；
- 重复3和4，直到$\Delta T<0.01 K$。
① ②两个过程并非简单的逆过程，否则2过程得到的$T_0’=T_0$，即迭代没有意义。实际上之所以这么做，是因为过程2在干绝热过程计算时，每一层混合气体的总$C_p$随着气体组分变化。并且，单一组分气体的$C_p$也在随温度变化的，而过程1使用了上下恒定的$C_p$。 下图展示了最初几步① ② ③过程示意图：
代码如下：
// estimate surface temperature and pressure
Real Ps = P0 * exp(-x1min / H0);
Real Ts = T0 * pow(Ps / P0, Rd / cp);
while (iter++ < max_iter) {
// read in vapors
air.w[iH2O] = xH2O;
air.w[iNH3] = xNH3;
air.w[IPR] = Ps;
air.w[IDN] = Ts;
// stop at just above P0
for (int i = is; i <= ie; ++i) {
pthermo->Extrapolate(&air, -dlnp / 2.,
Thermodynamics::Method::DryAdiabat);
if (air.w[IPR] < P0) break;
}
// extrapolate down to where air is
pthermo->Extrapolate(&air, log(P0 / air.w[IPR]),
Thermodynamics::Method::DryAdiabat);
// make up for the difference
Ts += T0 - air.w[IDN];
if (std::abs(T0 - air.w[IDN]) < 0.01) 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
air.ToMoleFraction();
for (int k = ks; k <= ke; ++k)
for (int j = js; j <= je; ++j) {
air.SetZero();
air.w[iH2O] = xH2O;
air.w[iNH3] = xNH3;
air.w[IPR] = Ps;
air.w[IDN] = Ts;
int i = is;
for (; i <= ie; ++i) {
// check relative humidity
Real rh = pthermo->RelativeHumidity(air, iNH3);
air.w[iNH3] *= std::min(rh_max_nh3 / rh, 1.);
AirParcelHelper::distribute_to_primitive(this, k, j, i, air);
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
if (air.w[IDN] < Tmin) break;
}
// Replace adiabatic atmosphere with isothermal atmosphere if temperature
// is too low
pthermo->Extrapolate(&air, dlnp, Thermodynamics::Method::DryAdiabat);
for (; i <= ie; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::Isothermal);
AirParcelHelper::distribute_to_primitive(this, k, j, i, air);
}
}
手动修正大气地不稳定度，使用adlnTdlnP，adlnNH3dlnP这个输入量来完成，会作用到Extrapolate函数上。
// if requested
// modify atmospheric stability
Real adlnTdlnP = pin->GetOrAddReal("problem", "adlnTdlnP", 0.);
Real adlnNH3dlnP = pin->GetOrAddReal("problem", "adlnNH3dlnP", 0.);
if (adlnTdlnP != 0.) {
Real pmin = pin->GetOrAddReal("problem", "adlnTdlnP.pmin", 1.);
Real pmax = pin->GetOrAddReal("problem", "adlnTdlnP.pmax", 20.);
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);
auto &&air = AirParcelHelper::gather_from_primitive(this, 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);
}
}
}
if (adlnNH3dlnP != 0.) {
Real pmin = pin->GetOrAddReal("problem", "adlnNH3dlnP.pmin", 1.);
Real pmax = pin->GetOrAddReal("problem", "adlnNH3dlnP.pmax", 20.);
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);
auto &&air = AirParcelHelper::gather_from_primitive(this, k, j, ibegin);
air.ToMoleFraction();
for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
air.w[iNH3] += adlnNH3dlnP * air.w[iNH3] * dlnp;
auto rates = pthermo->TryEquilibriumTP_VaporCloud(air, iNH3);
air.w[iNH3] += rates[0];
AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
}
}
}
设置痕量气体和电子，碱金属。计算蒸汽压和浓度：
// set tracers, electron and Na
int ielec = pindex->GetTracerId("e-");
int iNa = pindex->GetTracerId("Na");
auto ptracer = pimpl->ptracer;
for (int k = ks; k <= ke; ++k)
for (int j = js; j <= je; ++j)
for (int i = is; i <= ie; ++i) {
Real temp = pthermo->GetTemp(this, k, j, i);
Real pH2S = xH2S * phydro->w(IPR, k, j, i);
Real pNa = xNa * phydro->w(IPR, k, j, i);
Real svp = sat_vapor_p_Na_H2S_Visscher(temp, pH2S);
pNa = std::min(svp, pNa);
ptracer->u(iNa, k, j, i) = pNa / (Constants::kBoltz * temp);
ptracer->u(ielec, k, j, i) = saha_ionization_electron_density(
temp, ptracer->u(iNa, k, j, i), 5.14);
}
// primitive to conserved conversion (hydro)
peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
js, je, ks, ke);
// conserved to primitive conversion (tracer)
peos->PassiveScalarConservedToPrimitive(pscalars->s, phydro->u, pscalars->r,
pscalars->r, pcoord, is, ie, js, je,
ks, ke);
最后是辐射传输计算。
// Microwave radiative transfer needs temperatures at cell interfaces, which
// are interpolated from cell centered hydrodynamic variables. Normally, the
// boundary conditions are taken care of internally. But, since we call
// radiative tranfer directly in pgen, we would need to update the boundary
// conditions manually. The following lines of code updates the boundary
// conditions.
phydro->hbvar.SwapHydroQuantity(phydro->w, HydroBoundaryQuantity::prim);
pbval->ApplyPhysicalBoundaries(0., 0., pbval->bvars_main_int);
// calculate radiative transfer
auto prad = pimpl->prad;
for (int k = ks; k <= ke; ++k) {
// run RT models
app->Log("run microwave radiative transfer");
for (int j = js; j <= je; ++j) prad->CalRadiance(this, k, j);
}
本节涉及到了使用湿绝热垂直递减率来构建大气温度廓线的方法，考虑了多种水凝物潜热释放的情形，具体请参考文献【Moist Adiabats with Multiple Condensing Species: A New Theory with Application to Giant-Planet Atmospheres】