LBNL - Review, Fix, and Verify Simulator Calculations

[ecosang]

Below is the PR message using LaTeX equations, grounded in the key equations from _get_interior_cv_temp_estimate. Copy everything between the horizontal rules into your GitHub PR description (GitHub renders $...$ and $$...$$ LaTeX in Markdown).

---

Fix: Duplicate specific-heat factor in the thermal storage term of TFSimulator

Summary

The tensorized energy-balance update in tf_simulator.py multiplies the specific heat capacity $c$ twice in the thermal storage (energy-accumulation) term, producing a $c^2$ factor instead of $c$. This over-weights the storage term by a factor of $\sim c$ (≈ 1000 for building air), corrupting the transient temperature response. This PR removes the duplicate multiplication in both the numerator and the denominator.

Reference: the base interior-node equation

The scalar reference implementation in _get_interior_cv_temp_estimate documents the correct interior-CV energy balance. With uniform spacing ($u = v = \delta_x$) and uniform conductivity ($k_1 = k_2 = k_3 = k_4 = k$), the balance is:

$k(vz)\frac{T_{i-1,j}-T_{i,j}}{u}+ k(uz)\frac{T_{i,j-1}-T_{i,j}}{v}+ k(vz)\frac{T_{i+1,j}-T_{i,j}}{u}+ k(uz)\frac{T_{i,j+1}-T_{i,j}}{v}+ Q_x= \frac{\rho c, u v z}{\Delta t}\left(T_{i,j}-T_{i,j}^{(-)}\right)$

The right-hand side is the energy storage term. It derives from the time derivative of internal energy $E = m c T = \rho,(uvz),c,T$:

$$ \frac{dE}{dt} = \rho,(uvz),c,\frac{dT}{dt} ;\approx; \frac{\rho c, u v z}{\Delta t}\left(T_{i,j}-T_{i,j}^{(-)}\right) $$

The specific heat $c$ appears exactly once. A units check confirms it, and shows it is consistent with the conduction coefficient $\tfrac{k A_n}{\delta_x}$:

$$ \left[\frac{\rho c, uvz}{\Delta t}\right] = \frac{\tfrac{\text{kg}}{\text{m}^3}\cdot\tfrac{\text{J}}{\text{kg}\cdot\text{K}}\cdot \text{m}^3}{\text{s}} = \frac{\text{J/K}}{\text{s}} = \frac{\text{W}}{\text{K}} = \left[\frac{k A_n}{\delta_x}\right] $$

Solving for $T_{i,j}$ (matching the documented closed form) gives:

$$ T_{i,j} = \frac{\displaystyle\sum_{\text{neighbors}} T_{\text{neighbor}} + \frac{Q_x}{zk} + t_0, T_{i,j}^{(-)}}{4 + t_0}, \qquad t_0 = \frac{\rho c, \delta_x^2}{k,\Delta t} = \frac{\delta_x^2}{\Delta t,\alpha}, \quad \alpha = \frac{k}{\rho c} $$

Tensorized form and the storage coefficient

TFSimulator uses the volumetric (non-normalized) form of the same balance. Define the storage coefficient:

$$ S \equiv \frac{\rho c, U V z}{\Delta t} $$

so the tensor update is (schematically, with $N_0, D_0$ the conduction/convection/mass/$Q_{\text{lwx}}$ terms that do not contain $c$):

$$ T = \frac{N_0 + S, T^{(-)}}{D_0 + S} $$

The bug replaces $S$ with $c,S = \dfrac{\rho c^2 U V z}{\Delta t}$ in both numerator and denominator.

Why the $c^2$ error does not cancel

Because $N_0$ and $D_0$ do not contain $c$, the erroneous factor does not cancel:

$$ T_{\text{correct}} = \frac{N_0 + S,T^{(-)}}{D_0 + S}, \qquad T_{\text{buggy}} = \frac{N_0 + cS,T^{(-)}}{D_0 + cS} $$

For building air $c \approx 1000$, the storage term dominates and drives:

$$ \lim_{c\to\infty} T_{\text{buggy}} = \frac{cS,T^{(-)}}{cS} = T^{(-)} $$

i.e. each CV barely updates per iteration — the transient response is slowed by $\sim c$. (The converged steady state is still correct, since there $T = T^{(-)}$, but the dynamics and per-step behavior are wrong.)

The bug in code

Denominator — _get_denominator:

# BEFORE (buggy): c multiplied twice -> c^2
dt3 = tf.math.multiply(t_density, self._t_u)        # rho * U
dt3 = tf.math.multiply(dt3, self._t_v)              # rho * U * V
dt3 = tf.math.multiply(dt3, t_heat_capacity)        # rho * U * V * c
dt3 = tf.scalar_mul(t_z, dt3)                       # rho * U * V * c * z
dt3 = tf.math.multiply(dt3, t_heat_capacity)        # rho * U * V * c^2 * z  <-- BUG
dt3 = tf.math.divide(dt3, t_delta_t)

Numerator — _get_numerator:

# BEFORE (buggy): c multiplied twice -> c^2
nt3 = tf.math.multiply(t_density, self._t_u)        # rho * U
nt3 = tf.math.multiply(nt3, self._t_v)              # rho * U * V
nt3 = tf.math.multiply(nt3, t_heat_capacity)        # rho * U * V * c
nt3 = tf.scalar_mul(t_z, nt3)                       # rho * U * V * c * z
nt3 = tf.math.multiply(nt3, t_heat_capacity)        # rho * U * V * c^2 * z  <-- BUG
nt3 = tf.math.multiply(nt3, t_temp_minus)           # ... * T^(-)
nt3 = tf.math.divide(nt3, t_delta_t)

The fix

Remove the second tf.math.multiply(..., t_heat_capacity) in both functions so the storage term equals $\dfrac{\rho c, U V z}{\Delta t}$ (and $\times T^{(-)}$ in the numerator):

# AFTER (correct): rho * c * U * V * z / dt
dt3 = tf.math.multiply(t_density, self._t_u)
dt3 = tf.math.multiply(dt3, self._t_v)
dt3 = tf.math.multiply(dt3, t_heat_capacity)
dt3 = tf.scalar_mul(t_z, dt3)
dt3 = tf.math.divide(dt3, t_delta_t)

# AFTER (correct): rho * c * U * V * z / dt * T^(-)
nt3 = tf.math.multiply(t_density, self._t_u)
nt3 = tf.math.multiply(nt3, self._t_v)
nt3 = tf.math.multiply(nt3, t_heat_capacity)
nt3 = tf.scalar_mul(t_z, nt3)
nt3 = tf.math.multiply(nt3, t_temp_minus)
nt3 = tf.math.divide(nt3, t_delta_t)

Verification

• Corrected storage coefficient $\dfrac{\rho c, UVz}{\Delta t}$ has units W/K, matching the conduction/convection coefficients — dimensionally consistent.
• Matches the scalar reference balance in _get_interior_cv_temp_estimate and the docstring of update_temperature_estimates, which already stated the storage term as $\tfrac{C\rho UVz}{\Delta t}$ (the docstring was correct; only the code was wrong).
• File parses/compiles cleanly.

---

In Consultation With / Contributions By: Yun-Dam Ko @yun-dam

[ecosang]

@s2t2 based on discussion with Yun-dam @yun-dam

[ecosang]

@s2t2 I have this typo in the docstring. the code is fine. This is the only issue. But I updated the docstring a bit more comprehensive.

[s2t2]

ok thanks!

@ecosang or @yun-dam can you please post a screenshot of the updated docs site content?

@amcberkes let us know if the formula updates look good to you. I believe some related changes may have been recently made internally?

[ecosang]

this is the revised one. mainly $q_{lwx} $ is W/m2 but miss to multiply area.

**this is previous. **

[amcberkes]

@s2t2 yes we have made that change internally in March.

[ecosang]

@StS2 seems several synchronization issues here. How do you want to deal with this? I may close this PR if Google's repo already has this. During that time, I will work with @yun-dam to find the other potential issues. Any comment regarding the syncing process? do you want me to merge this PR #169 to open-source one? but at least, I want to modify the "docs" update. I can update it separately.

[s2t2]

@ecosang I would be inclined to merge this PR. We will include it in the content that needs to get synced, and process it in the same way.

Will take a final look on Monday and probably merge if everything looks good.