CenterForTheBuiltEnvironment / pythermalcomfort

Package to calculate several thermal comfort indices (e.g. PMV, PPD, SET, adaptive) and convert physical variables.
https://pythermalcomfort.readthedocs.io/en/latest/
MIT License
147 stars 52 forks source link

The differences in two_nodes between R and Python #143

Open TwinGan opened 3 weeks ago

TwinGan commented 3 weeks ago

Differences caused by additional code logic unique to Python

  1. First At the end of the loop, Py package will reassign the two variables e_skin, m_rsw( regsw in R) When w > w_max, the e_rsw will be assign to prsw * emax which will make m_rsw become different because the reassign logic.
    
        if w > w_max:
            w = w_max
            p_rsw = w_max / 0.94
            e_rsw = p_rsw * e_max
            e_diff = 0.06 * (1.0 - p_rsw) * e_max
        if e_max < 0:
            e_diff = 0
            e_rsw = 0
            w = w_max
        e_skin = (
            e_rsw + e_diff
        )  # total evaporative heat loss sweating and vapor diffusion
        m_rsw = (
            e_rsw / 0.68
        )  # back calculating the mass of regulatory sweating as a function of e_rsw

In the R package, there was no this logic.

[The extra logic](https://github.com/CenterForTheBuiltEnvironment/pythermalcomfort/blob/496f3799de287737f2ea53cc6a8c900052a29aaa/pythermalcomfort/models/two_nodes.py#L421C8-L426C85)

2. Second
For input:
calc2Node(ta = 45,tr = 45,vel = 1.1,rh = 20,clo = 0.2,met = 3) in R
two_nodes(tdb = 45,tr = 45,v = 1.1,rh = 20,clo = 0.2,met = 3) in Python

The R code result will lead Python's result by one loop
It 's the result of m_rsw(regsw) in the end of each loop (I have added the following extra logic to R)

![image](https://github.com/user-attachments/assets/818bd78a-c77a-4ffb-a34d-6b733eb52e89)
![image](https://github.com/user-attachments/assets/05c49c17-62dd-4c03-8fbd-b7548045b950)

You can see that the 59 loop m_rsw in R equals the 60 loop in Python.
I am sure that these two loop all start from 1
![image](https://github.com/user-attachments/assets/7747817b-5eb1-42f3-b42d-770620763db8)
![image](https://github.com/user-attachments/assets/aa8f3ee9-f832-41d4-a25a-15d2e1df408b)

I don't understand why but this situation happened.

3. Could we increase the tolerance of W. 
Expected w = 0.2 Actual w = 0.226183272947463
Expected w = 0.2 Actual w = 0.174338146533957
Expected w = 0.1 Actual w = 0.127748749502237
All three cases failed because of the tolerance value.

@FedericoTartarini and @marcelschweiker  what shall I do with those differences?
marcelschweiker commented 3 weeks ago

@TwinGan and @FedericoTartarini If I interprete correctly, we can fix the difference either by making an addition to R or removing the addition in py. I am fine with both and suggest we go for the version, that is closer to the reference values. @FedericoTartarini What was the source again for the reference values for 2-node model?

FedericoTartarini commented 2 weeks ago

@TwinGan and @marcelschweiker please find below my comments.

First point

The code from Gagge is

image

The two-node code from Mark Fountain is as follows.

image

I have also attached both documents for your reference. Please let me know what do you think.

Gagge et al_1986_A standard predictive index of human reponse to the thermal environemnt.pdf

Fountain_A derivation of the Gagge 2-node model.pdf

Second point

I am not sure why this is happening. It must be something related to how we initialise the simulation. One must be skipping the first run.

Third point

I think there is an issue in pythermalcomfort and we need to return w with more decimals, I think I have already partially fixed this in the new version of pythermalcomfort. I should not be returning w with only one dwecimal value. Please have a look at the code in development. If the code is not there yet I will fix it as soon as I have time.

FedericoTartarini commented 1 week ago

@TwinGan, Marcel and I looked into the code and we have noticed that the code is actually the same. The difference you mentioned in point 1 is not there. However, we identified an issue caused by the iteration and this is what it is causing the error.

I will fix point 2 and 3 in the new version of pythermalcomfort.