AllInOne for ISC25 Coding Challenge

并非笔记. 实际上是流水账.

Reading

三个主要循环? 读一下.

Data Structures & Parameters

请问您今天要来点数据吗?

in graupel.cpp and graupel.hpp

graupel() 的输入列表.

/**
 * @brief
 *
 * @param [in] nvec Number of horizontal points
 * @param [in] ke Number of grid points in vertical direction
 * @param [in] ivstart Start index for horizontal direction
 * @param [in] ivend End index for horizontal direction
 * @param [in] kstart Start index for vertical direction
 * @param [in] dt Time step for integration of microphysics (s)
 * @param [in] dz Layer thickness of full levels (m)
 * @param [inout] t Temperature in Kelvin
 * @param [in] rho Density of moist air (kg/m3)
 * @param [in] p Pressure (Pa)
 * @param [inout] qv Specific water vapor content (kg/kg)
 * @param [inout] qc Specific cloud water content (kg/kg)
 * @param [inout] qi Specific cloud ice content (kg/kg)
 * @param [inout] qr Specific rain content (kg/kg)
 * @param [inout] qs Specific snow content  kg/kg)
 * @param [inout] qg Specific graupel content (kg/kg)
 * @param [in] qnc Cloud number concentration
 * @param [out] prr_gsp Precipitation rate of rain, grid-scale (kg/(m2*s))
 * @param [out] pri_gsp Precipitation rate of ice, grid-scale (kg/(m2*s))
 * @param [out] prs_gsp Precipitation rate of snow, grid-scale (kg/(m2*s))
 * @param [out] prg_gsp Precipitation rate of graupel, grid-scale (kg/(m2*s))
 * @param [out] pflx Total precipitation flux
 *
 */

参数名	类型	描述	单位
nvec	[in]	水平网格点的数量	-
ke	[in]	垂直方向网格点的数量 (高度)	-
ivstart	[in]	水平方向起始索引	-
ivend	[in]	水平方向结束索引 (结束是 0)	-
kstart	[in]	垂直方向起始索引	-
dt	[in]	微物理过程积分的时间步长	秒 (s)
dz	[in]	完整层的层厚	米 (m)
t	[inout]	温度	开尔文 (K)
rho	[in]	湿空气的密度	kg/m³
p	[in]	气压	帕 (Pa)
qv	[inout]	比水汽含量	kg/kg
qc	[inout]	比云水含量	kg/kg
qi	[inout]	比云冰含量	kg/kg
qr	[inout]	比雨水含量	kg/kg
qs	[inout]	比雪含量	kg/kg
qg	[inout]	比冰雹含量	kg/kg
qnc	[in]	云中粒子数浓度	-
prr_gsp	[out]	网格尺度降雨率	kg/(m²·s)
pri_gsp	[out]	网格尺度冰降率	kg/(m²·s)
prs_gsp	[out]	网格尺度降雪率	kg/(m²·s)
prg_gsp	[out]	网格尺度冰雹降率	kg/(m²·s)
pflx	[out]	总降水通量	-

说明:

• [in] 表示只读取参数。
• [inout] 表示既读取又修改参数。
• [out] 表示仅修改参数。
• 表格中列出了每个参数的名称、类型（输入/输出属性）、描述以及单位（若适用）。

t_qx_ptr 和 q 分别存储了 6 种形态的水的比含量, 其中 4 种形态的水需要预测降水概率 (凝结水)?

struct t_qx_ptr {
  t_qx_ptr(array_1d_t<real_t> &p_, array_1d_t<real_t> &x_) : p(p_), x(x_) {} 
  array_1d_t<real_t> &p;
  array_1d_t<real_t> &x;
}; 

  q.emplace_back(prr_gsp, qr);
  q.emplace_back(pri_gsp, qi);
  q.emplace_back(prs_gsp, qs);
  q.emplace_back(prg_gsp, qg);
  q.emplace_back(emptyArray, qc);
  q.emplace_back(emptyArray, qv);

第一个元素 p 是降<水>率, 是 2D 的, 长度为水平方向的地块数量 nvec, 也就是说只预测地面是否降水; 第二个元素 x 是比水凝物质量, 是 3D 的, 长度为空间内所有地块的数量nvec × ke, 也就是说有不同高度的区别.

初始:

id	p	x	explain
0	prr_gsp	qr	雨
1	pri_gsp	qi	冰
2	prs_gsp	qs	雪
3	prg_gsp	qg	冰雹
4	emptyArray	qc	云
5	emptyArray	qv	水汽

可以这样记忆:

• pr<预测水种类>_gsp
  • Precipitation Rate, 沉淀率
• q<水种类>
  • Quantity, 量

还有一些定义在 graupel() 中定义的本地变量/向量.

  array_1d_t<bool> is_sig_present(nvec *
                                  ke); // is snow, ice or graupel present?

  array_1d_t<size_t> ind_k(nvec * ke),
      ind_i(nvec * ke); // k index of gathered point, iv index of gathered point
  array_2d_t<size_t> kmin(
      nvec, array_1d_t<size_t>(np)); // first level with condensate

  real_t cv, vc, eta, zeta, qvsi, qice, qliq, qtot, dvsw, dvsw0, dvsi, n_ice,
      m_ice, x_ice, n_snow, l_snow, ice_dep, e_int, stot, xrho;

  real_t update[3], // scratch array with output from precipitation step
      sink[nx],     // tendencies
      dqdt[nx];     // tendencies
  array_1d_t<real_t> eflx(
      nvec); // internal energy flux from precipitation (W/m2 )
  array_2d_t<real_t> sx2x(nx, array_1d_t<real_t>(nx, ZERO)), // conversion rates
      vt(nvec,
         array_1d_t<real_t>(
             np)); // terminal velocity for different hydrometeor categories

  // 中间省略 q 的定义, 此事在上文已有记载

  size_t jmx = 0;
  size_t jmx_ = jmx;

向量, 其中意义列由注释中总结而来, 一维数组的 * 表示乘法运算, 而二维数组的 × 用于分割第一维和第二维的长度.

变量名	类型	大小	用途/意义	可能单位	用于	循环 (含线程安全)
is_sig_present	一维布尔	nvec * ke	表示是否存在雪、冰或霰	无（布尔值）	1	下标依赖 jmx
ind_k	一维整数	nvec * ke	收集点的垂直索引	无（索引）	1	下标依赖 jmx
ind_i	一维整数	nvec * ke	收集点的水平索引	无（索引）	1	下标依赖 jmx
kmin	二维整数	nvec × np	存在凝结物的第一个垂直层级	无（索引）	1	✅ 需最值规约
update[3]	一维实数	3	降水步骤的输出暂存数组	视上下文而定	3	✅ 可私有
sink[nx]	一维实数	nx	倾向量，可能表示水汽或质量的损失率	kg/(kg·s)	2	✅ 可私有
dqdt[nx]	一维实数	nx	倾向量，可能是水汽混合比的变化率	kg/(kg·s)	2	✅ 可私有
eflx	一维实数	nvec	降水引起的内部能量通量	W/m²	2, 3	✅ eflx[iv] 可私有
sx2x	二维实数	nx × nx	转换率，如液态到固态的水相变化率	s⁻¹	2	✅ 线程私有
vt	二维实数	nvec × np	不同水凝物类别的末端速度	m/s	3	✅ vt[iv] 可私有

其中的标量值 (用途推测):

变量名	类型/结构	用途/意义	可能单位
cv	real_t	可能表示某种比热容	J/(kg·K)
vc	real_t	可能表示某种速度或体积相关量	m/s 或 m³
eta	real_t	可能是动力粘度或效率因子	kg/(m·s) 或 无
zeta	real_t	可能是某种阻尼系数或中间变量	无
qvsi	real_t	冰面上的饱和比湿	kg/kg
qice	real_t	冰的质量混合比	kg/kg
qliq	real_t	液态水的质量混合比	kg/kg
qtot	real_t	总水汽或水物质的质量混合比	kg/kg
dvsw	real_t	水面上的水汽扩散系数	m²/s
dvsw0	real_t	水面上的参考水汽扩散系数	m²/s
dvsi	real_t	冰面上的水汽扩散系数	m²/s
n_ice	real_t	冰粒子的数浓度	m⁻³
m_ice	real_t	冰的质量	kg
x_ice	real_t	冰的某种特征尺寸或质量相关变量	m 或 kg
n_snow	real_t	雪粒子的数浓度	m⁻³
l_snow	real_t	雪的某种长度或质量相关变量	m 或 kg
ice_dep	real_t	冰的沉积率	kg/(m²·s)
e_int	real_t	内部能量	J/kg
stot	real_t	总熵或某种总量	J/(kg·K)
xrho	real_t	可能是密度相关变量	kg/m³
jmx	size_t	符合条件格点计数器	无

除了 jmx 外都是线程私有的, 而 jmx 若需要可进行求和规约.

in common.hpp

namespace idx {
constexpr size_t nx = 6;  // number of water species
constexpr size_t np = 4;  // number of precipitating water species
constexpr size_t lqr = 0; // index for rain
constexpr size_t lqi = 1; // index for ice
constexpr size_t lqs = 2; // index for snow
constexpr size_t lqg = 3; // index for graupel
constexpr size_t lqc = 4; // index for cloud
constexpr size_t lqv = 5; // index for vapor

constexpr size_t qx_ind[] = {lqv, lqc, lqr, lqs, lqi, lqg};
constexpr size_t qp_ind[] = {lqr, lqi, lqs, lqg};
} // namespace idx

在 common.hpp 中有 namespace idx, 定义了水种类数量和下标及顺序, 预测水种类数量和下标及顺序.

Functions

我的函数很多.

in core/properti/common.hpp

/**
 * @brief Calculates saturation pressure over water
 *
 * @param [in] t Temperature (kelvin)
 * @return Saturation pressure
 */
TARGET real_t sat_pres_water(real_t t) {
  return c1es * exp(c3les * (t - thermodyn::tmelt) / (t - c4les));
}

/**
 * @brief Saturation vapor pressure (over ice) at constant density
 *
 * @param [in] t Temperature (kelvin)
 * @param [in] rho Density
 * @return saturation pressure
 */
TARGET real_t qsat_ice_rho(real_t t, real_t rho) {
  return sat_pres_ice(t) / (rho * rv * t);
}

$$ e_{\mathrm{sat}}(T) = c_{1es} \cdot \exp\left( \frac{c_{3les} \cdot (T - T_{\text{melt}})}{T - c_{4les}} \right) \ q_{\mathrm{sat}}(T, \rho) = \frac{e_{\mathrm{sat}}(T)}{\rho \cdot R_v \cdot T} $$

给定温度 t 和气体密度 rho, 计算冰面上的饱和比湿, 即单位质量空气中水汽的最大含量. 而冰面上饱和蒸汽比湿 $q_{\mathrm{sat}}$ 依赖于液态水面上的饱和蒸气压 $e_{\mathrm{sat}}$.

1st Loop

我的循环很少.

  // The loop is intentionally i<nlev; since we are using an unsigned integer
  // data type, when i reaches 0, and you try to decrement further, (to -1), it
  // wraps to the maximum value representable by size_t.
  size_t oned_vec_index;
  for (size_t i = ke - 1; i < ke; --i) {
    for (size_t j = ivstart; j < ivend; j++) {
      oned_vec_index = i * ivend + j;
      if ((std::max({q[lqc].x[oned_vec_index], q[lqr].x[oned_vec_index],
                     q[lqs].x[oned_vec_index], q[lqi].x[oned_vec_index],
                     q[lqg].x[oned_vec_index]}) > qmin) or
          ((t[oned_vec_index] < tfrz_het2) and
           (q[lqv].x[oned_vec_index] >
            qsat_ice_rho(t[oned_vec_index], rho[oned_vec_index])))) {
        jmx_ = jmx_ + 1;
        ind_k[jmx] = i;
        ind_i[jmx] = j;
        is_sig_present[jmx] =
            std::max({q[lqs].x[oned_vec_index], q[lqi].x[oned_vec_index],
                      q[lqg].x[oned_vec_index]}) > qmin;
        jmx = jmx_;
      }

      for (size_t ix = 0; ix < np; ix++) {
        if (i == (ke - 1)) {
          kmin[j][qp_ind[ix]] = ke + 1;
          q[qp_ind[ix]].p[j] = ZERO;
          vt[j][ix] = ZERO;
        }

        if (q[qp_ind[ix]].x[oned_vec_index] > qmin) {
          kmin[j][qp_ind[ix]] = i;
        }
      }
    }
  }

第一个 for 从 ke - 1 到 0 枚举, 第二个 for 从 ivstart 到 ivend 枚举格点编号.

内层:

• 根据"某种条件"决定是否处理, 记录 ind_k, ind_i, is_sig_president.
  • "某种条件":
    • 若任意一种形态的水的比含量 q[lq<水种类>].x[<格点下标>] > qmin;
    • 或着:
      • 格点温度 t[<格点下标>] < tfrz_het2;
      • 且水含量 q[lqv].x[<格点编号>] > 冰面上饱和比热湿 qsat_ice_rho(t[<格点编号>], rho[<格点编号>]).
• 记录 kmin[j][qp_ind[ix]], 其中 qp_ind[ix] 是要预测的水.
  • 若是顶点, 即 i == (ke - 1):
    • 则初始概率 q[qp_ind[ix]].p[j] 和某种速度 vt[j][ix] 置 0;
  • 并且记录每种凝结物质存在的最低的层级 kmin[j][qp_ind[ix]] = i.

2nd Loop

我独自循环.

水凝物相变, 质量守恒过程.

3rd Loop

末日三环.

计算结果, 模拟沉降, 并且计算因沉降导致的能量输出.

Thinking

我思故我在.

简单观察发现, 似乎局部变量和输入的向量都是线程安全的, 至少在 nvec 的角度上来看.

那么理解后面两个循环具体做了什么似乎就不是特别重要, 这可以放在后面具体优化的时候来做. 我们只需要知道这个程序进行了三个操作:

• 标记需要计算的点
• 计算这些点
• 更新每个格点的另一些东西.

更进一步分析, 如果不在乎第二个循环的负载均衡的话, 可以每个线程分治水平方向的地块 [ivstart, ivend), 并负责全部垂直高度 ke.

Coding

男子程序员的日常.

Analyzing

理科生尝试沉入分析.

串行代码计时

串行亭奇谭

单位是 ms.

数据集/循环	1st	2nd	3rd
11k	1	3	3
20k	204	1162	783

Journal

Space: the final frontier. These are the voyages of the starship Enterprise. Its five-year mission: to explore strange new worlds; to seek out new life and new civilizations; to boldly go where no man has gone before!

Dynamic To-Do List

空.

Day 3: 04/04

阅读代码, 和生佳诺开会确认现在进度.