Implementation of Monte Carlo method

This page details the implementation of Monte Carlo method, see [Modest] and [Ozturk] for more details.

The code are stored in project : github.com/feelpp/hygrothermie

In folder src/rht (github.com/feelpp/hygrothermie/tree/master/src/rhtt)

1. Algorithm in pseudo-code

For this work, I choose three algorithms (the principle difference is the order of execution of different element).

1.1. Element

This section present different function that use the global algorithm.

1.1.1. Generate Source

\begin{algorithm}
    \caption{Source Point}
    \begin{algorithmic}
    \FUNCTION{IsInTriangle}{$\mathcal{T}$, $Q$}
    \STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$
    \STATE $u \leftarrow  \frac{(\mathbf{e_2} \cdot \mathbf{e_2}) (\vec{PQ} \cdot \mathbf{e_1}) - (\mathbf{e_1} \cdot \mathbf{e_2}) (\vec{PQ} \cdot \mathbf{e_2}) }{(\mathbf{e_1} \cdot \mathbf{e_1})(\mathbf{e_2} \cdot \mathbf{e_2})-(\mathbf{e_1} \cdot \mathbf{e_2})^2}$
    \STATE $v \leftarrow \frac{(\mathbf{e_1} \cdot \mathbf{e_1}) (\vec{PQ} \cdot \mathbf{e_2}) - (\mathbf{e_1} \cdot \mathbf{e_2}) (\vec{PQ} \cdot \mathbf{e_1}) }{(\mathbf{e_1} \cdot \mathbf{e_1})(\mathbf{e_2} \cdot \mathbf{e_2})-(\mathbf{e_1} \cdot \mathbf{e_2})^2}$
    \RETURN $(0 \leq u \leq 1)(0 \leq v \leq 1)$
    \ENDFUNCTION

    \FUNCTION{SourcePoint}{$\mathcal{T}$}
    \STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$

    \STATE $\xi_1$ and $\xi_2$ random number in $[0,1]$
    \STATE $\phi \leftarrow 2\pi \xi_1$
    \STATE $r \leftarrow R \sqrt{\xi_2}$
    \STATE $Q \leftarrow \begin{pmatrix} x_C + r cos \phi \\ y_C + r sin \phi \end{pmatrix}$

    \WHILE{ IsInTriangle($\mathcal{T}$, $Q$)}
        \STATE $\xi_1$ and $\xi_2$ random number in $[0,1]$
        \STATE $\phi \leftarrow 2\pi \xi_1$
        \STATE $r \leftarrow R \sqrt{\xi_2}$
        \STATE $Q \leftarrow \begin{pmatrix} x_C + r cos \phi \\ y_C + r sin \phi \end{pmatrix}$
    \ENDWHILE
    \RETURN $Q$
    \ENDFUNCTION
    \end{algorithmic}
\end{algorithm}

1.1.2. Generate Ray

\begin{algorithm}
    \caption{GenerateRay}
    \begin{algorithmic}
    \FUNCTION{Ray}{$\mathcal{T}$, $s$}
    \STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$
    \STATE $s$ source point

    \STATE $\phi_n \leftarrow $ incident angle of $normal(\mathcal{T})$
    \STATE $\theta_n \leftarrow $ azimuthal angle of $normal(\mathcal{T})$

    \STATE $\phi_{random} \leftarrow$ random number in $[0,\pi]$
    \STATE $\theta_{random} \leftarrow$ random number in $[0,2\pi]$

    \STATE $\phi \leftarrow \phi_n + \phi_{random}$
    \STATE $\theta \leftarrow \theta_n + \theta_{random}$

    \RETURN $\phi$, $\theta$
    \ENDFUNCTION
    \end{algorithmic}
\end{algorithm}

1.1.3. Visibility

\begin{algorithm}
    \caption{Visibility}
    \begin{algorithmic}
    \FUNCTION{Visibility}{$\mathcal{T}_1$, $\mathcal{T}_2$}

    \STATE $\mathbf{n_1} \leftarrow$ $normal(\mathcal{T}_1)$
    \STATE $\mathbf{n_2} \leftarrow$ $normal(\mathcal{T}_2)$

    \RETURN $\left( \mathbf{n_1} \cdot \mathbf{n_2} < 0 \right)$
    \ENDFUNCTION
    \end{algoritmic}
\end{algorithm}

1.1.4. Interception

\begin{algorithm}
    \caption{Interception}
    \begin{algorithmic}

    \FUNCTION{Interception}{$r$,$\mathcal{T}$}

    \STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$
    \STATE $r$ is ray, compose by origin point $O$ and directional vector $\mathbf{v}$

    \STATE $div = \left( \mathbf{v} \wedge \mathbf{e_2} \right) \cdot \mathbf{e_1} \ne 0$

    \IF{$div \ne 0$}
        \STATE $u = \frac{ \vec{OP}  \cdot \left( \mathbf{e_2} \wedge \mathbf{v} \right)}{div}$
        \IF{$0 \leq u \leq 1$}
            \STATE $t = \frac{ \vec{OP}  \cdot \left( \mathbf{e_1} \wedge \mathbf{v} \right)}{div}$
            \IF{$0 \leq t \leq 1$}
                \RETURN 1
            \ENDIF
        \ENDIF
    \ENDIF

    \RETURN 0
    \ENDFUNCTION

    \end{algorithmic}
\end{algorithm}

1.2. First Strategy

Here, for all element \(e^i_k\) of a surface \(A_i\), we generate once a set of rays and for all element \(e^i_l\), we browse the set of ray to know which ray is intercepted.

\begin{algorithm}
    \caption{ViewFactor}
    \begin{algorithmic}
    \FUNCTION{ViewFactor1}{$A_i, A_j, N_s, N_r$}
    \STATE $A_i$ is surface decomposed by $K$ triangles $e^i_k$
    \STATE $A_j$ is surface decomposed by $L$ triangles $e^j_l$
    \STATE $N_s$ number of source point of one element
    \STATE $N_r$ number of ray from one source point

    \STATE $Area_i \leftarrow 0$ : area of surface $A_i$
    \STATE $\Sigma_2 \leftarrow 0$

    \FOR{$e^i_k$ in set of element of surface $A_i$}

        \FOR{$1 \leq i_s \leq N_s$}
            \STATE $s \leftarrow$ \textbf{SourcePoint}($\mathcal{T}$)
            \FOR{$1 \leq i_r \leq N_r$}
                \STATE $r \leftarrow$ \textbf{GenerateRay}($\mathcal{T}$, $r$)
            \ENDFOR
        \ENDFOR

        \STATE $\Sigma_1 \leftarrow 0$
        \FOR{$e^j_l$ in set of element of surface $A_j$}
            \IF{\textbf{Visibility}($e^i_k$,$e^j_l$)}
                \STATE $n \leftarrow 0$ : number of ray intercept
                \FOR{ray $r$ in set of ray}
                    \IF{\textbf{Interception}($r$,$e^j_k$)}
                        \STATE $n \leftarrow n + 1$
                    \ENDIF
                \ENDFOR
                \STATE $F_{local} \leftarrow \frac{n}{N_s N_r}$
                \STATE $\Sigma_1 \leftarrow \Sigma_1 + F_{local}$
            \ENDIF
        \ENDFOR
        \STATE $Area_i \leftarrow Area_i + Area(e^i_k)$
        \STATE $\Sigma_2 \leftarrow \Sigma_2 + Area(e^i_k) \Sigma_1$
    \ENDFOR
    \STATE $F_{global} \leftarrow \frac{1}{Area(A_i)} \Sigma_2$
    \RETURN $F_{global}$
    \ENDFUNCTION
    \end{algorithmic}
\end{algorithm}

Complexity : \(O\left( K (L+1) N_s N_r \right) = O\left( K L N_s N_r \right)\), with \(K\) number of elements of \(A_i\), \(L\) number of elements of \(A_j\), \(N_s\) number of source points and \(N_r\) number of rays from one source point.

Complexity more precise :

\begin{eqnarray} K \times \left(N_s \times complexity(\mathbf{SourcePoint}) \times N_r \times complexity(\mathbf{Ray}) \right. \\ \left. + L \times complexity(\mathbf{Visibility}) \times N_r \times N_s \times complexity(\mathbf{Interception}) \right) \end{eqnarray}

Disadvantage : requires more memories. (not really important)

1.3. Second Strategy

Here, for all element \(e^i_k\) of surface \(A_i\), we generate once a set of rays and for all ray, we browse the set of element \(e^j_l\) to know which element intercept the ray.

\begin{algorithm}
    \caption{ViewFactor}
    \begin{algorithmic}
    \FUNCTION{ViewFactor2}{$A_i, A_j, N_s, N_r$}
    \STATE $A_i$ is surface decomposed by $K$ triangles $e^i_k$
    \STATE $A_j$ is surface decomposed by $L$ triangles $e^j_l$
    \STATE $N_s$ number of source point of one element
    \STATE $N_r$ number of ray from one source point

    \STATE $Area_i \leftarrow 0$ : area of surface $A_i$
    \STATE $\Sigma_2 \leftarrow 0$

    \FOR{$e^i_k$ in set of element of surface $A_i$}
        \STATE $n \leftarrow \{n_l\}_{1\leq l \leq L}$ is a list of zeros
        \FOR{$1 \leq i_s \leq N_s$}
            \STATE $s \leftarrow$ \textbf{SourcePoint}($\mathcal{T}$)
            \FOR{$1 \leq i_r \leq N_r$}
                \STATE $r \leftarrow$ \textbf{GenerateRay}($\mathcal{T}$, $s$)
                \STATE $l \leftarrow 0$
                \WHILE{\textbf{Visibility}($e^i_k$,$e^j_l$) and not\textbf{Interception}($r$,$e^j_k$) and $l \leq L$}
                    \STATE $l \leftarrow l+1$
                \ENDWHILE
                \IF{$l \leq L$}
                    \STATE $n_j \leftarrow n_j+1$
                \ENDIF
            \ENDFOR
        \ENDFOR
        \STATE $Area_i \leftarrow Area_i + Area(e^i_k)$
        \STATE ($F_{e^i_k \rightarrow e^j_l} \leftarrow  \frac{n_j}{N_s N_r}$)
        \STATE $F_{e^i_k \rightarrow A_j} \leftarrow \frac{\Sigma(n_j)}{N_s N_r} \mbox{ or } = \Sigma(F_{e^i_k \rightarrow e^j_l} )$
        \STATE $\Sigma_2 \leftarrow \Sigma_2 + Area(e^i_k) F_{local} $
    \ENDFOR
    \STATE $F_{global} \leftarrow F_{global} + \frac{1}{Area(A_i)} \Sigma_2$
    \RETURN $F_{global}$
    \ENDFUNCTION
    \end{algorithmic}
\end{algorithm}

Complexity : Worst case \(O\left( KN_sN_r \frac{L(L+1)}{1} \right)\), with \(K\) number of elements of \(A_i\), \(L\) number of elements of \(A_j\), \(N_s\) number of source points and \(N_r\) number of rays from one source point.

1.4. Third Strategy

Here, for each couple of elements \(e^i_l\), \(e^j_k\), we calculate the corresponding view factor. This algorithm comes from Ozturk.

\begin{algorithm}
    \caption{Calculation of View Factor}
    \begin{algorithmic}
    \FUNCTION{ViewFactor3}{$A_i$, $A_j$, $N_s$, $N_r$}
    \STATE $A_i$ is surface decomposed by $K$ triangles $e^i_k$
    \STATE $A_j$ is surface decomposed by $L$ triangles $e^j_l$
    \STATE $N_s$ number of source point of one element
    \STATE $N_r$ number of ray from one source point

    \STATE $Area_i \leftarrow 0$ : area of surface $A_i$
    \STATE $\Sigma_2 \leftarrow 0$

    \FOR{$e^i_k$ in set of element of surface $A_i$}
        \STATE $\Sigma_1 \leftarrow 0$
        \FOR{$e^j_l$ in set of element of surface $A_j$}
            \IF{\textbf{Visibility}($e^i_k$,$e^j_l$)}
                \STATE $n \leftarrow 0$
                \FOR{$1 \leq i_s \leq N_s$}
                    \STATE $s \leftarrow$ \textbf{SourcePoint}($\mathcal{T}$)
                    \FOR{$1 \leq i_r \leq N_r$}
                        \STATE $r \leftarrow$ \textbf{GenerateRay}($\mathcal{T}$,$s$)
                        \IF{\textbf{Interception}($r$,$e^j_k$)}
                            \STATE $n \leftarrow n + 1$
                        \ENDIF
                    \ENDFOR
                \ENDFOR
                \STATE $F_{e^i_k \rightarrow e^j_l} \leftarrow \frac{n}{N_s N_r}$
            \ENDIF
            \STATE $\Sigma_1 \leftarrow \Sigma_1 + F_{e^i_k \rightarrow e^j_l}$
        \ENDFOR
        \STATE $Area_i \leftarrow Area_i + Area(e^i_k)$
        \STATE $\Sigma_2 \leftarrow \Sigma_2 + Area(e^i_k) \Sigma_1$
    \ENDFOR
    \STATE $F_{global} \leftarrow \frac{1}{Area(A_i)} \Sigma_2$
    \RETURN  $F_{global}$
    \ENDFUNCTION
    \end{algorithmic}
\end{algorithm}

Complexity : \(O\left( K L N_s N_r \right)\), with \(K\) number of elements of \(A_i\), \(L\) number of elements of \(A_j\), \(N_s\) number of source points and \(N_r\) number of rays from one source point.

Complexity more precise :

\[\begin{eqnarray} K \times L \times N_s \times complexity(\mathbf{SourcePoint}) \times N_r \\ \times complexity(\mathbf{Ray}) \times complexity(\mathbf{Interception})] \end{eqnarray}\]

1.5. Choose the algorithm

The choose depends on the complexity of SourcePoint.

2. Environment of development

The application radiative-heat-transfert (which calculates the heat transfert and viex factor) is developed in software Feelpp. Feelpp is a software of modelisation (heat, solid mechanics, fluid mechanics …) [feelpp].

3. Structutre of code

The source codes are place in src/rht :

rht_option.cpp : the option of radiative-heat-transfert application
rht_main.cpp : the main script of radiative-heat-transfert application

I haven’t time to begining the codes with Feel++. The main structure (rht_main.cpp) is writed by M. Prudhomme.

The folder rough_draft regroups the codes not terminated.

References

[Incropera] Fundamentals of Heat and Mass Transfer, Theodore L. Bergman, Adrienne S. Lavine, Frank P. Incropera, David P. DeWitt, Wiley.
[Kumaran1996] Kumaran, M. Kumar, Heat, Air and Moisture Transfer in Insulated Envelope Parts, Final Report, Volume 3.
[Ozturk] Abdurrahman Ozturk, Implementation of View Factor Model and Radiative Heat Transfer Model in MOOSE. University of South Carolina, Download PDF
[Modest] M.F. Modest, Radiative Heat Transfet. Elsevier Science, 2013, ISBN 9780123869906, Chapter 8
[Visionray] Visionaray: A Cross-Platform Ray Tracing Template Library, Zellmann, Stefan and Wickeroth, Daniel and Lang, Ulrich, Proceedings of the 10th Workshop on Software Engineering and Architectures for Realtime Interactive Systems (IEEE SEARIS 2017), 2017
[Visionary::github] Visionaray, Zellman, github.com/szellmann/visionaray
Méthode Numérique pour le EDP, Christophe Prudhomme, feelpp.github.io/csmi-edp/#/, Programmation de le méthode élément fini, Formule de quadrature
CEMOSIS, CEMOSIS, www.cemosis.fr/
L’Accord de Paris, unfccc.int/fr/process-and-meetings/the-paris-agreement/l-accord-de-paris, UNFCC Sites and platforms
[Infographie] Énergie : quel secteur d’activité consomme le plus ?, www.gazprom-energy.fr/gazmagazine/2017/09/consommation-energetique-secteur-activite/, Gazprom Energy
L’électricité dans le secteur résidentiel, www.edf.fr/groupe-edf/espaces-dedies/l-energie-de-a-a-z/tout-sur-l-energie/le-developpement-durable/l-electricite-dans-le-secteur-residentiel, EDF
Synapse Concept, www.synapse-concept.com/
4fastsim-ibat, www.cemosis.fr/projects/4fastsim-ibat/
Feel++ Docs, docs.feelpp.org/feelpp/index.html