# Projet de 4fastsim (method of PBDW)

## 2. PBDW’s method

### 2.1. Introduction of PBDW’s method

Abstract: The principal objective of this part is to study a numerical tool which is less expensive for simulate the thermal distribution of the building. The construction a model thermique which

Here we introduce a non-intrusive reduced order method of data assimilation for parameterized PDEs: the PBDW (Parameterized-Background Data-Weak). These method belong to the family of Reduced Basis methods.

## 3. Methodologie

### 3.1. Stationary problem of the generic form

$\mathcal{P} : \Omega \times \mathcal{D} \rightarrow \mathbb{R}$

where $\mathcal{P}$ is the problem,$\Omega \in \mathbb{R}^d$ is the physical domain of dim $d = 2,3$ , $\mathcal{D} \in \mathbb{R}^{Np}$ is the parameter domain.

We will consider the solution $u \in \mathcal{X}$, a Hilbert space s.t $H^1_0(\Omega) \in \mathcal{X} \in H^1(\Omega)$ .

We assume we have some knowledge of the system encoded in the parameterized PDE model $\mathcal{P}^{bk}$ (i.e the best adapted model available for the problem $\mathcal{P})$, resulting in the best knowlege state $u^{bk}(\textbf{p}) \in \mathcal{X}$.

### 3.2. Sensor

We also assume we have a second source of information in the form of field observations.

We have $M$ experimental observations, we assume our data $y^{obs}_m, 1 \leq m \leq M$, we desire output quantity

$y_m^{obs}(\textbf{p}) = l_m(u^{true}(\textbf{p}))$

where $u^{true}(\textbf{p})$ represents the true physical state of the system and $l_m$ are linear functionals representing the sensors

### 3.3. Goal

$\Rightarrow$ The challenge is fomulate a scheme to interpret and employ data from mersurements in some optimal way to contribute to the knowledge of the state in the numerical model

$\Rightarrow$ In order to combine the best knowledge of the system modeled by $\mathcal{P}^{bk}$ and the experimental observations of the physical state, we will consider a PBDW solution of two parts: background and update.

$\Rightarrow$ The goal is to fomulate the PBDW method in which we will approximate the true physical state $u^{true}(\textbf{p})$ by

$$$u^{bk}(\textbf{p}) + \eta$$$

where $\eta \in \mathcal{X}$ is an update correction term associated to the experimental observations.

### 3.4. Reduced basis background

The idea of reduced basis methods is to compute an inexpensive and accurate approximation, $u^{bk}_N(\textbf{p})$ of $\mathcal{P}^{bk}$ for any $\textbf{p} \in \mathcal{D}^{bk}$ by seeking a linear combination of particular solutions

$$$u^{bk}_N(\textbf{p}) = \sum^{N}_{i=1}\beta_i(\textbf{p}) u^{bk}_N(\textbf{p}_i)$$$

We perform a Gram-Schmidt orthnormalization, and introduce the Background space $\mathcal{Z}_N = span\{\zeta_i\}^{N}_{i=1} \subset \mathcal{X}$

### 3.5. Data-informed update

We will consider the experimental observations to contribute in the update piece.

Use the the representation

$\varphi_m = exp\left( \frac{-(x - x_m)^2}{2r^2} \right) s.t \int_{\omega} \varphi_m(x) d \Omega =1$

where $x_m \in \mathbb{R}^d$ is the center of the $m^{th}$ sensor of radius $r$.

To evaluate the information these sensors can gather from a physical state $v \in \mathcal{X}$, we define the function

$l_m(v) = \int_{\Omega} \varphi_m(x) v d \Omega \quad 1 \leq m \leq M$

We need to use (3) to construct the second approximation space $\mathcal{U}_M \subset\mathcal{X}$ in which we find the update $\eta$.

Let us define the Riesz operator $\mathcal{R}_{\mathcal{X}}: \mathcal{X}' \rightarrow \mathcal{X}$ s.t

$<v,\mathcal{R}_{\mathcal{X}} l_m>_{\mathcal{X}} = l(v)$

Introduce the update basis functions $q_m = \mathcal{R}_{\mathcal{X}} l_m$ and for any physical state $u^{true}$ of the configuration, we have

$<u^{true},\mathcal{R}_{\mathcal{X}} q_m>_{\mathcal{X}} = l_m(u^{true})$

We define Update Space $\mathcal{U}_{M_{max}} = span \{q_m \}^{M_{max}}_{m=1}$ where $M_{max}$ is the maximum number of sensors available.

### 3.6. PBDW problem statement

The PBDW aims at approximating the true physical state utrue(p) for some configuration $\textbf{p}$ by

$u_{N,M} = z_M + \eta_M$

where $z_N$ is in $\mathcal{Z}_N$, corresponds RB approximation pf best knowledge solution $u^{bk}(\textbf{p})$

and the term $\eta_M$ is in $\mathcal{U}_M$, is the correction term associated with $M$ observations

We pose the PBDW approximation as the solution to the following minimization problem. Find $(u_{N,M} \in \mathcal{X}, z_N \in \mathcal{Z}_N,\eta_M \in \mathcal{U}_M)$ s.t

$(U_{N,M};Z_{N};\eta_{M})=\underset{U_{N,M}\in \mathcal{X},Z_{N}\in \mathcal{Z}_{N},\eta_{M}\in \mathcal{U}_{M}}{arginf}\{ \| \eta_{M}\|_{\mathcal{X}}^{2}\} \\ <U_{N,M}-Z_{N},v>_{\mathcal{X}}=<\eta_{M},v>_{X}\ \forall\ v\in \mathcal{X} \\ <U_{N,M},\phi>_{\mathcal{X}}=<u(p),\phi>_{\mathcal{X}} \forall\ \phi \in \mathcal{U}_{M}$

Find $(\eta_M()\textbf{p} \in \mathcal{U}_M, z_N(\textbf{p})\in \mathcal{Z}_N)$ s.t

$<\eta_M,q>_{\mathcal{X}} + <z_N,q>_{\mathcal{X}} = <u^{true}(\textbf{p},q)>_{\mathcal{X}} \quad q \in \mathcal{U}_M$
$<\eta_M,p>_{\mathcal{X}} = 0 \quad \forall p \in \mathcal{Z}_N$

Use the definition of update basis functions $q_m \in \mathcal{X}$ we have

$<\eta_M,q_m>_{\mathcal{X}} = y^{obs}_m(\textbf{p})$
$y^{obs}_m(\textbf{p}) = l_m(u^{true}(\textbf{p}))$

The corresponding algebraic formulation is to find $\eta_M \in \mathbb{R}^M, z_N \in \mathbb{R}^N$ s.t

$\begin{pmatrix} A&B\\ B^{T}&0\\ \end{pmatrix} \begin{pmatrix} \overrightarrow{\eta_{M}}\\ \overrightarrow{Z_{N}}\\ \end{pmatrix} =\begin{pmatrix} \overrightarrow{y^{obs}}\\ 0\\ \end{pmatrix}$

where $(y^{obs}_m = y^{obs}_m)$ , $\textbf{A}_{m,m} = <q_m,q_m>_{\mathcal{X}}$ and $\textbf{B}_{m,n} = <\zeta_n, \zeta_m>_{\mathcal{X}}$

For a more stable numerical implementation, we can consider solving :

$\textbf{B}^T \textbf{A}^{-1} z_N = \textbf{B}^T \textbf{A}^{-1} y^{obs}$
$\eta_M = \textbf{A}^{-1} (y^{obs} - \textbf{B} z_N)$

*Offline stage:

1. Choose $N$, $M$.

2. Choose $\mathcal{P}^{bk}$ and $\mathcal{D}^{bk}$ .

3. Construct $\mathcal{Z}_N$ and $\mathcal{U}_M$.

4. Construct $\textbf{A}$, $\textbf{B}$

*Online

1. Solve (3.17) and (3.18)

2. Calcule output

$l_m(u_N,u_M) = \sum^{M}_{m=1} n_M l_m(q_m) + \sum^{N}_{n=1} z_N l_m(\zeta_n)$

### 3.7. Algorithms

Parameters

Given : Sample

$Q = [Q_1, \cdots , Q_2]$

$T_{ext} = [T_1, \cdots , T_2]$

$h = [h_1, \cdots , h_2]$

$\sum_{test} = (\textbf{p}_1, \cdots, \textbf{p}_{ntest}) \in \mathcal{D}^{ntest}$

$SN = (textbf{p}_1 , \cdots, textbf{p}_N)$

Buil Background space $\mathcal{Z}_N [ \mathcal{N} \times N ]$

(Sans Greedy)

` runheat ` ($\textbf{p}_1$) $\mapsto u_h(\textbf{p}_1)$

$\mathcal{Z}_1 = span(u_h(\textbf{p}_1))$

for $\textbf{p}_i$ in $SN$ do:

$\quad \quad$ ` runheat `($\textbf{p}_i$) $\mapsto u_h(\textbf{p}_i)$

$\quad \quad$ $z_i = u_h(\textbf{p}_i) - \sum_{u \in \mathcal{Z}_{i-1}} proj_u (u_h(\textbf{p}_i))$ (where $proj_u (v) = \frac{<u,v>}{<u,u>}u$ )

$\quad \quad$ $\mathcal{Z}_i = \mathcal{Z}_{i-1} \cup \{\mathcal{Z}_i\}$

Normaliser: for $i = 1$ to $N$ do:

$\quad \quad$ $\mathcal{Z}_i = \frac{z_1}{<z_i, z_i>}$

(Avec Greedy)

Choose ramdomly $\textbf{p}_1 \in \mathcal{D}$

$S_1 = \{ \textbf{p}_1\}$ and $\mathcal{Z}_1 = span(u_h(\textbf{p}_1))$

for $i = 2$ to $N$ do:

$\quad \quad$ $\textbf{p}_i = \underset{\textbf{p} \in \sum_{test}}{argmax} \frac{|| u_h(/textbf{p}) - \sum_{u \in \mathcal{Z}_{i-1}} proj_u (u_h(\textbf{p})) ||_{L^2}}{ || u_h(\textbf{p}) ||_{L^2} }$

$\quad \quad$ $S_i = S_{i-1} \cup \{ \textbf{p}_i \}$

$\quad \quad$ $\mathcal{Z}_i = \mathcal{Z}_{i-1} \cup \{ u_h(\textbf{p}_i) \}$

Normaliser: for $i = 1$ to $N$ do:

$\quad \quad$ $\mathcal{Z}_i = \frac{z_1}{<z_i, z_i>}$

Build Update space $\mathcal{U}_M [\mathcal{N} \times M ]$

Define $l_m(v) \forall v$ test function in $\mathcal{X}_h$

for $m = 1$ to $M$ do:

$\quad \quad$ Define sensor functionals over ball $B((x_m,y_m), r)$

$\quad \quad$ $\varphi_m = exp \left( - \frac{(x - x_m)^2 (y - y_m)^2}{2r^2} \right) s.t \int_{w} \varphi_m (x) d \Omega = 1$

$\quad \quad$ Solve for $q_m : <v,q_m>_{H^1} = l_m(v)$

$\quad \quad$ $a(q_m,v) = \int_{\Omega} \nabla q_m \cdot \nabla v + \int_{\Omega} q_m v$

$\quad \quad$ $l_m(v) = \int_{\Omega} \varphi_m v$

$\quad \quad$ $\mathcal{U}_M [:,m ] = q_m$

Build matrix PBDW

$(\mathbb{A}_{i,j}) = <q_j,q_i>_{H^1} = a(q_j, q_i) \in \mathbb{R}^{M \times M}$

$(\mathbb{B}_{i,j}) = <\zeta_j,q_i>_{H^1} = a(\zeta_j, q_i) \in \mathbb{R}^{M \times N}$

$\mathbb{K}^{PBDW} =\begin{pmatrix} A&B\\ B^{T}&0\\ \end{pmatrix}$

Build $y^{obs}$

Choose randomly $\textbf{p} \in \mathcal{D}$

` runfluid ` $(\textbf{p}) \mapsto u^{true}(\textbf{p})$

for $m = 1$ to $M$:

$\quad \quad$ $y^{obs}_m(\textbf{p}) = l_m(u^{true}(\textbf{p}))$

$y^{obs}[:, m] = y^{obs}_m(\textbf{p})$

Solve PBDW

$\begin{pmatrix} \overrightarrow{\eta_{M}}\\ \overrightarrow{z_{N}}\\ \end{pmatrix} = (\mathbb{K}^{PBDW})^{-1} * \begin{pmatrix} \overrightarrow{y^{obs}}\\ 0\\ \end{pmatrix}$

Output

$u_{N,M} = \sum^{M}_{i=1} \eta_M^i q_i + \sum^{N}_{j=1} z^j_N \zeta_j$

### 3.8. Performance

In this case, the output quantity over the basis functions of the two approximation spaces can be precalculated, allowing for evaluation of the output of the PBDW approximation in $\mathcal{O}(N + M)$ operations, without fully reconstructing the PBDW approximation from the basis functions ${\zeta_n}^N_{n=1}$ and ${q_m}^M_m=1$, a procedure in $\mathcal{O}(\mathcal{N}_h)$ operations.

### 3.9. PBDW error and stability considerations

The question is how can we know if our model is well-posed ? The problem dependend obviously in the construction of the Background and Ipdate spaces. So we want to compute a value which can depend on the two approximation spaces. We introduce the inf-sup Stability constant .

Stability constant
$\beta_{N,M} = inf_{w \in \mathcal{Z}_N} sup_{v \in \mathcal{U}_M} \frac{<w,v>_{\mathcal{X}}}{ \| w \|_{\mathcal{X}} \| v \|_{\mathcal{X}}}$

Note that $\beta_{N,M}$ is an non-increasing function of $N$ and non-decreasing function of $M$ with $\beta_{N,M} = O$ for $N>M$.

Present the best-fit of the approximation spaces

Best-fit error
$\frac{inf_{q \in \mathcal{U}_M} inf_{z \in \mathcal{Z}_N} \| u^{true} - z - q \|_{\mathcal{X}}}{\| u^{true} \|_{\mathcal{X}}}$

A priori error estimatation is derived for the formulation as a function of the stability constant and the best-fit of the approximation spaces

Priory error
$\frac{\| u^{true} - u_{N,M} \|_{X}}{\| u^{true} \|_{\mathcal{X}}} \leq \left( 1 + \frac{1}{\beta_{N,M}} \right) \frac{ inf_{q \in \mathcal{U}_M} inf_{z \in \mathcal{Z}_N} \| u^{true} - z - q \|_{\mathcal{X}}}{\| u^{true} \|_{\mathcal{X}}}$

Note that this error estimate of PBDW method dependend strongly on the stability constant. So if we have obtion of choosing $M$, we need to build the approxiamation spaces in order to maximize the stability of the formular and also minimize the best-fit error.

On the other hand, in pratique, we mostly don’t have a true physical state $u^{true}$ to precise the priory error, so we work in the boundary of error such as the second term in the formular of PBDW error. Better approximation, smaller the boundary of error.

## 4. Numerical implementation of PBDW

We have to link 2 toolboxes of `Heat transfert` and `heatfluid` to solve the problem. For each toolboxes, we need to creat a `.json` file to precise the materials properties

Main.cpp
``````/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4*/

#include <feel/feelmodels/heatfluid/heatfluid.hpp>
#include <feel/feelmodels/heat/heat.hpp>
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <stdlib.h>
#include <complex>
#include <algorithm>
#include <cstdlib>

namespace Feel
{
static const int nDim = 2;
static const int OrderT = 1;
static const int OrderV = 2;
static const int OrderP = 1;
typedef FeelModels::Heat< Simplex<nDim,1>,
Lagrange<OrderT, Scalar,Continuous,PointSetFekete> > model_heat_type;

typedef FeelModels::FluidMechanics< Simplex<nDim,1>,
Lagrange<OrderV, Vectorial,Continuous,PointSetFekete>,
Lagrange<OrderP, Scalar,Continuous,PointSetFekete> > model_fluid_type;
typedef FeelModels::HeatFluid< model_heat_type,model_fluid_type> model_heatfluid_type;

template <typename HeatFluidType>
std::map<std::string,double>
runApplicationHeatFluid( boost::shared_ptr<HeatFluidType> heatFluid, std::map<std::string,double> const& params )
{
for ( auto const& param : params ) //param in params

heatFluid->solve();
heatFluid->exportResults();
if ( heatFluid->isStationary() )
{
heatFluid->solve();
heatFluid->exportResults();
}
else
{
#if 1
auto heatModel = heatFluid->heatModel();
auto fluidModel = heatFluid->fluidModel();
if ( true )
{
auto myexpr = expr( soption(_prefix=heatModel->prefix(),_name="initial-solution.temperature"),
"",heatModel->worldComm(),heatModel->repository().expr() );
heatModel->fieldTemperaturePtr()->on(_range=heatModel->rangeMeshElements(),_expr=myexpr);
}
fluidModel->fieldVelocityPressurePtr()->setConstant(0);
heatModel->timeStepBdfTemperature()->start( heatModel->fieldTemperature() );
fluidModel->timeStepBDF()->start( *fluidModel->fieldVelocityPressurePtr() );
heatModel->updateTime( heatFluid->heatModel()->timeStepBase()->time() );
fluidModel->updateTime( heatFluid->heatModel()->timeStepBase()->time() );
heatFluid->updateTime( heatFluid->heatModel()->timeStepBase()->time() );
fluidModel->blockVectorSolution().updateVectorFromSubVectors();
heatModel->blockVectorSolution().updateVectorFromSubVectors();
#endif
for ( ; !heatFluid->timeStepBase()->isFinished(); heatFluid->updateTimeStep() )
{
if (heatFluid->worldComm().isMasterRank())
{
std::cout << "============================================================\n";
std::cout << "time simulation: " << heatFluid->time() << "s \n";
std::cout << "============================================================\n";
}

heatFluid->solve();
//heatFluid->exportResults();
}
heatFluid->heatModel()->exportMeasures( 0 );
}

std::map<std::string,double> res;
for ( std::string const& s : std::vector<std::string>({ "temperature_pointA", "temperature_pointB" ,"temperature_pointC","temperature_pointD","temperature_pointE","temperature_pointF","temperature_pointG","temperature_pointH","temperature_pointI","temperature_pointK","temperature_pointL","temperature_pointM"}) )
res[s] = heatFluid->heatModel()->postProcessMeasuresIO().measure( s );
return res;
}

template <typename HeatType>
auto runApplicationHeat( boost::shared_ptr<HeatType> heat, std::map<std::string,double> const& params )
{
for ( auto const& param : params )

heat->solve();
heat->exportResults();
auto const& e = heat->fieldTemperature();

return e; //neu muon them dk kiem tra if thi can khai bao e o ngoai, vay no la type gi?
}

template <typename HeatType>
double produitscalair(boost::shared_ptr<HeatType> heat, model_heat_type::element_temperature_type const& u, model_heat_type::element_temperature_type const& v)
{
double Lg = 2.;
double p = integrate(_range=elements( heat->mesh()), _expr= inner(idv(u),idv(v))).evaluate()(0,0);
double val = g + Lg*p;
//double val = inner_product( u,v);
return val;
}

template <typename HeatType>
double produit(boost::shared_ptr<HeatType> heat, model_heat_type::element_temperature_type const& u, model_heat_type::element_temperature_type const& v)
{
double val = inner_product( u,v);
return val;
}

template <typename HeatType>
auto proj(boost::shared_ptr<HeatType> heat,model_heat_type::element_temperature_type const& u, model_heat_type::element_temperature_type const& v)
{
double uv = produit(heat, u, v);
double uu = produit(heat,u,u);
double scalar = uv / uu;
std::cout << "scalar = "<< scalar << "\n";
auto proju = heat->spaceTemperature()->element( scalar*idv(u) );
return proju;
}

} // namespace Feel

int myrandom (int i) { return std::rand()%i;}

struct Coordinate
{
double x, y;
}pointA,pointB,pointC,pointD,pointE,pointF,pointG,pointH,pointI,pointK,pointL,pointM;

int main( int argc, char** argv )
{
using namespace Feel;
int N;
int M;
double r;

po::options_description pbdwdatafromaerooptions( "application pbdw datafromaero options" );
("case.dimension", Feel::po::value<int>()->default_value( 3 ), "dimension")
("case.discretization", Feel::po::value<std::string>()->default_value( "P1-P2P1" ), "discretization : P1-P2P1")
("M", Feel::po::value<int>(&M)->default_value( 1), "nbObservation")
("N", Feel::po::value<int>(&N)->default_value( 2), "nbParameters")
("r", Feel::po::value<double>(&r)->default_value( 0.15), "rayon")
;

Environment env( _argc=argc, _argv=argv,
_desc=pbdwdatafromaerooptions,
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));

boost::shared_ptr<model_heatfluid_type> heatFluid( new model_heatfluid_type("heat-fluid") );
heatFluid->init();
heatFluid->printAndSaveInfo();

boost::shared_ptr<model_heat_type> heat( new model_heat_type("heat",false) );
heat->setMesh( heatFluid->mesh() );
heat->setStationary( true );
heat->init();
heat->printAndSaveInfo();

tic();

//Creat Sample SN
std::vector<std::map<std::string,double>> SN;
std::map<std::string,double> params0;
params0["T_heater"] = 310;
params0["T_ext"] = 280;
params0["h"] = 1.0/(0.06+0.01/0.5 + 0.3/0.8 + 0.20/0.032 +0.016/0.313 +0.14);
SN.push_back(params0 );

std::map<std::string,double> params1;
params1["T_heater"] = 300;
params1["T_ext"] = 280;
params1["h"] = 1.0/(0.06+0.01/0.5 + 0.3/0.8 + 0.20/0.032 +0.016/0.313 +0.14);
SN.push_back(params1);

std::map<std::string,double> params2;
params2["T_heater"] = 300; //330
params2["T_ext"] = 290; //310
params2["h"] = 0.1;
SN.push_back(params2);

//std::vector<std::map<std::string,double>> SN;
//for(int i=0; i< N; i++)
//  SN[i] = DN[i];

using element_temperature_type = typename model_heat_type::element_temperature_type;
std::vector< element_temperature_type > ZN;

for (int i=0; i<N; i++)
{
auto params = SN[i];
auto u = runApplicationHeat( heat, params );
ZN.push_back(u);
//std::cout << "Zn " << ZN.back().size() << " vs " << u.size() << "\n";
//std::cout << "Zn " << ZN.back().sum() << " vs " << u.sum() << "\n";
}

//Build Background space ZN_orthn/////////////////////////////////////////////////////////////////////////
std::vector< element_temperature_type > ZN_orth( ZN.size(),heat->spaceTemperature()->element() );
ZN_orth[0] = ZN[0];
for(int i = 1; i < N; i++)
{
auto s = heat->spaceTemperature()->element();
for(int j = 0; j < i; j++)
{
s += proj(heat, ZN_orth[j], ZN[i]);
}
ZN_orth[i] = ZN[i] - s;
}
for(int i = 0; i < ZN_orth.size(); i++)
{
double norm = std::sqrt(produitscalair(heat, ZN_orth[i], ZN_orth[i]));
ZN_orth[i].scale( 1./norm);
}

//check
for (int i = 0; i < ZN_orth.size(); i++)
{
for (int j = 0; j < ZN_orth.size(); j++)
{
double scar = produitscalair(heat, ZN_orth[i], ZN_orth[j]);
if ( Environment::isMasterRank() )
std::cout << "i="<<i<<"; j="<<j<<"; scar=" << scar << "\n";
}
}

std::vector< element_temperature_type > ZN_orthn( ZN_orth.size(),heat->spaceTemperature()->element() );
ZN_orthn[0] = ZN_orth[0];
for(int i = 1; i < N; i++)
{
auto s = heat->spaceTemperature()->element();
for(int j = 0; j < i; j++)
{
s += proj(heat, ZN_orthn[j], ZN_orth[i]);
}
ZN_orthn[i] = ZN_orth[i] - s;
}
for(int i = 0; i < ZN_orthn.size(); i++)
{
double norm = std::sqrt(produitscalair(heat, ZN_orthn[i], ZN_orthn[i]));
ZN_orthn[i].scale( 1./norm);
}

//check ZnOrthn?
for (int i = 0; i < ZN_orthn.size(); i++)
{
for (int j = 0; j < ZN_orthn.size(); j++)
{
double scar = produitscalair(heat, ZN_orthn[i], ZN_orthn[j]);
if ( Environment::isMasterRank() )
std::cout << "i="<<i<<"; j="<<j<<"; scar=" << scar << "\n";
}
}

//Buil update space UM////////////////////////////////////////////////////////////////////////////////////////
//Creat List of points
std::vector<Coordinate> points;
pointA.x = 0.05;
pointA.y = 0.25;
points.push_back(pointA);

pointB.x = 0.55;
pointB.y = 0.25;
points.push_back(pointB);

pointC.x = 1.05;
pointC.y = 0.25;
points.push_back(pointC);

pointD.x = 1.55;
pointD.y = 0.25;
points.push_back(pointD);

pointE.x = 0.05;
pointE.y = 0.5;
points.push_back(pointE);

pointF.x = 0.55;
pointF.y = 0.5;
points.push_back(pointF);

pointG.x = 1.05;
pointG.y = 0.5;
points.push_back(pointG);

pointH.x = 1.55;
pointH.y = 0.5;
points.push_back(pointH);

pointI.x = 0.05;
pointI.y = 0.75;
points.push_back(pointI);

pointK.x = 0.55;
pointK.y = 0.75;
points.push_back(pointK);

pointL.x = 1.05;
pointL.y = 0.75;
points.push_back(pointL);

pointM.x = 1.55;
pointM.y = 0.75;
points.push_back(pointM);

std::vector<int> data;
for (int i=0; i<M; ++i)
data.push_back(i);

std::random_shuffle ( data.begin(), data.end(), myrandom);

std::cout << "data contains:";
for (std::vector<int>::iterator it=data.begin(); it!=data.end(); ++it)
std::cout << ' ' << *it;
std::cout << '\n';

int p = rand() % N;
std::cout << " p = "<< p << "\n";
auto resHeatFluid = runApplicationHeatFluid( heatFluid, SN[p] );
auto u_true = heatFluid->heatModel()->fieldTemperature();
std::vector< double > y_obs;

auto phi = heat->spaceTemperature()->element( );
auto capteur = heat->spaceTemperature()->element( );

std::vector< element_temperature_type > UM( M,heat->spaceTemperature()->element() );
for(int i=0; i< M; i++)
{
int index = data[i];
auto sensor = points[index];
std::cout << "xm="<< sensor.x<< "; ym="<<sensor.y<< "\n";
auto phiExpr = exp( - (pow(Px()-sensor.x,2) + pow(Py()-sensor.y,2))/(2*r*r) )*chi( (pow(Px()-sensor.x,2) + pow(Py()-sensor.y,2))<r*r );

phi.on( _range=elements( heat->mesh()),_expr=phiExpr);
capteur += phi;

double c = integrate(_range=elements( heat->mesh()), _expr= phiExpr).evaluate()(0,0);
std::cout << " Is that c = 1 ? " << c << "\n";

double y = integrate(_range=elements( heat->mesh()), _expr= phiExpr * idv(u_true)).evaluate()(0,0);
y_obs.push_back(y);

auto qm = heat->spaceTemperature()->element("qm" );
auto v = heat->spaceTemperature()->element( );

auto l = form1( _test= heat->spaceTemperature());
l = integrate(_range=elements( heat->mesh()), _expr= phiExpr * id(v));

auto a = form2( _trial=heat->spaceTemperature(), _test=heat->spaceTemperature());
a += integrate(_range=elements( heat->mesh()), _expr= inner(idt(qm),id(v)));

a.solve(_rhs=l,_solution=qm);
UM[i] = qm;
std::cout << " UM " << UM[i].size() << " vs " << qm.size() << "\n";
std::cout << " UM " << UM[i].sum() << " vs " << qm.sum() << "\n";
}

//Build matrice PBDW//////////////////////////////////////////////////////////////////////////
//int M = coords.size();
//int N = ZN_orthn.size();
using Eigen::MatrixXd;
MatrixXd A(M,M);
for(int i=0; i<M; i++)
{
for(int j=0; j<M; j++)
A(i,j) = produitscalair(heat,UM[j], UM[i]);
}
std::cout << "A = \n" << A << "\n";

MatrixXd B(M,N);
for(int i=0; i<M; i++)
{
for(int j=0; j<N; j++)
B(i,j) = produitscalair(heat,ZN_orthn[j], UM[i]);
}
std::cout << "B = \n"  << B << "\n";

MatrixXd K(M+N,M+N);
K.setZero();
K.block(0,0,M,M) = A;
K.block(0,M,M,N) = B;
K.block(M,0,N,M) = B.transpose();
std::cout << "K = \n"  << K << "\n";

//Solve PBDW//////////////////////////////////////////////////////////////////////////
MatrixXd Y(M+N,1);
Y.setZero();
for(int i=0; i<M; i++)
Y(i,0)=y_obs[i];
std::cout << "Y = \n"  << Y << "\n";

toc("offline");

Eigen::EigenSolver<MatrixXd> esK(K);
std::cout << "The eigenvalues of K are:" << endl << esK.eigenvalues() << endl;
std::cout << "The matrix of eigenvectors, V, is:" << endl << esK.eigenvectors() << endl << endl;

std::vector< double > lambdaK;
for(int i=0; i< N+M; i++   )
{
lambdaK.push_back(esK.eigenvalues()[i].real());
std::cout << "lambda_K = \n"  << lambdaK.back() << "\n";
}

double lambdaK_min = *std::min_element(lambdaK.begin(), lambdaK.end());
double lambdaK_max = *std::max_element(lambdaK.begin(), lambdaK.end());
double cond_K = std::abs(lambdaK_max)/std::abs(lambdaK_min);
std::cout << "Cond(K) = "  << cond_K << "\n";

tic();
MatrixXd S(M+N,1);
S = K.inverse()*Y;
std::cout << "S = \n"  << S << "\n";

//Output//////////////////////////////////////////////////////////////////////////
MatrixXd eta(M,1);
//eta = S.block(0,0,M,0);
for (int i=0; i<M; i++)
eta(i,0)=S(i,0);
MatrixXd z(N,1);
//z = S.block(M,0,N,0);
for (int i=0; i<N; i++)
z(i,0)=S(i+M,0);

std::cout << "z = \n"  << z << "\n";
std::cout << "eta = \n"  << eta << "\n";

auto u_MN = heat->spaceTemperature()->element();
for(int i=0; i<M; i++)
u_MN += heat->spaceTemperature()->element(eta(i,0)*idv(UM[i]));

for(int i=0; i<N; i++)
u_MN += heat->spaceTemperature()->element(z(i,0)*idv(ZN_orthn[i]));

toc("online");
////Stability coefficient//////////////////////////////////////////
tic();
MatrixXd C(N,N);
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
C(i,j) = produitscalair(heat,ZN_orthn[i], ZN_orthn[j]);
}
std::cout << "C = \n"  << C << "\n";

MatrixXd K1(N,N);
K1 = (B.transpose() * A.inverse() * B ) * C.inverse();
std::cout << "K1 = \n"  << K1 << "\n";

Eigen::EigenSolver<MatrixXd> es(K1);
std::cout << "The eigenvalues of K1 are:" << endl << es.eigenvalues() << endl;
std::cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;

std::vector< double > lambda;
for(int i=0; i< N; i++   )
{
lambda.push_back(es.eigenvalues()[i].real());
//std::cout << "lambda = \n"  << lambda.back() << "\n";
}

double lambda_min = *std::min_element(lambda.begin(), lambda.end());
double beta = std::sqrt(lambda_min);
//std::cout << "lambda_min = "  << lambda_min << "\n";
std::cout << "beta = "  << beta << "\n";

/////priori error/////////////////////////////////////////////////////////
double l2 = normL2(_range=elements(heat->mesh()), _expr=idv(u_true)-idv(u_MN) );
//std::cout << "h1 = "  << h1 << "\n";
double c_l2 = normL2(_range=elements(heat->mesh()), _expr=idv(u_true) );
//std::cout << "c_l2 = "  << c_l2 << "\n";
double e_l2 = (1./c_l2) * l2;
std::cout << "err_l2 = "  << e_l2 << "\n";
//std::cout << "c_h1 = "  << c_h1 << "\n";
double e_h1 = (1./c_h1) * h1;
std::cout << "err_h1 = "  << e_h1 << "\n";

///boundary error///////////////////////////////////////////////////////
std::vector< double > err_h1;
std::vector< double > err_l2;
for(int i=0; i<N; i++)
{
for(int j=0; j<M; j++)
{
err_h1.push_back(normH1(_range=elements(heat->mesh()),

err_l2.push_back(normL2(_range=elements(heat->mesh()), _expr=idv(u_true)-idv(ZN_orthn[i]) - idv(UM[j]) ));
}
}
double min_h1 = *std::min_element(err_h1.begin(), err_h1.end());
double min_l2 = *std::min_element(err_l2.begin(), err_l2.end());

double best_fit_l2 = (1./c_l2) * min_l2;
std::cout << "best_fit_l2 = "  << best_fit_l2 << "\n";
double bound_l2 = (1 + 1./beta) * best_fit_l2;
std::cout << "bound_l2 = "  << bound_l2 << "\n";

double best_fit_h1 = (1./c_h1) * min_h1;
std::cout << "best_fit_h1 = "  << best_fit_h1 << "\n";
double bound_h1 = (1 + 1./beta) * best_fit_h1;
std::cout << "bound_h1 = "  << bound_h1 << "\n";

toc("beta");
///save solution///////////////////////////////////////////////////////////////
auto g = exporter( _mesh=heat->mesh(),_name="u_true" );
g->save();

auto u_err = heat->spaceTemperature()->element();
for(int i = 0; i<u_true.size(); i++)
u_err[i] = std::log(std::abs(u_true[i] - u_MN[i]));
auto er = exporter( _mesh=heat->mesh(), _name="u_err");
er->save();

auto d = exporter( _mesh=heat->mesh(),_name="u_h" );
for(int i = 0; i < ZN.size(); i++)
d->save();

auto q = exporter( _mesh=heat->mesh(),_name="qm" );
for(int i = 0; i < UM.size(); i++)
q->save();

auto n = exporter( _mesh=heat->mesh(),_name="Znorthn" );
for(int i = 0; i < ZN_orthn.size(); i++)
n->save();

auto e = exporter( _mesh=heat->mesh(), _name="u_MN");
e->save();

auto QM = heat->spaceTemperature()->element();
for(int i=0; i<UM.size(); i++)
QM += UM[i];
auto f = exporter( _mesh=heat->mesh(),_name="qm_total" );
f->save();

auto gg = exporter( _mesh=heat->mesh(),_name="sensor_location" );
gg->save();

auto u_velocity = heatFluid->fluidModel()->fieldVelocity();
auto h = exporter( _mesh=heat->mesh(),_name="u_velocity" );
h->save();

std::cout << "============================================================\n";

return 0;
}``````
aero.json
``````// -*- mode: javascript -*-
{
"Name": "Thermo dynamics",
"ShortName":"ThermoDyn",
"Parameters":
{
"T_heater":310,
"T_ext":280,
"h":0.1
},
"Models":
{
"use-model-name":1,
"fluid":
{
"equations": "Navier-Stokes"
}
},
"Materials":
{
"air":
{
"markers":["air"],
"physics":["fluid","heat"],
"rho":"1",
"mu":"2.65e-2",
"k":"0.03",
"Cp":"1004",
"beta":"0.003660" //0.00006900
}
},
"BoundaryConditions":
{
"velocity":
{
"Dirichlet":
{
"exterior-walls": { "expr":"{0,0}" },
"heater1": { "expr":"{0,0}" }
}
},
"temperature":
{
"Dirichlet":
{
"heater1": { "expr":"T_heater:T_heater" }
},
"Robin":
{
"exterior-walls":
{
"expr1":"h:h",// h coeff
"expr2":"T_ext:T_ext"// temperature exterior
}
}
}
},
"PostProcess":
{
"use-model-name":1,
"heat-fluid":
{
"Exports":
{
"fields":["fluid.velocity","fluid.pressure","heat.temperature","fluid.pid"]
}
},
"fluid":
{
},
"heat":
{
"Measures":
{
"Points":
{
"pointA":
{
"coord":"{0.05,0.25}",
"fields":"temperature"
},
"pointB":
{
"coord":"{0.55,0.25}",
"fields":"temperature"
},
"pointC":
{
"coord":"{1.05,0.25}",
"fields":"temperature"
},
"pointD":
{
"coord":"{1.55,0.25}",
"fields":"temperature"
},
"pointE":
{
"coord":"{0.05,0.5}",
"fields":"temperature"
},
"pointF":
{
"coord":"{0.55,0.5}",
"fields":"temperature"
},
"pointG":
{
"coord":"{1.05,0.5}",
"fields":"temperature"
},
"pointH":
{
"coord":"{1.55,0.5}",
"fields":"temperature"
},
"pointI":
{
"coord":"{0.05,0.75}",
"fields":"temperature"
},
"pointK":
{
"coord":"{0.55,0.75}",
"fields":"temperature"
},
"pointL":
{
"coord":"{1.05,0.75}",
"fields":"temperature"
},
"pointM":
{
"coord":"{1.55,0.75}",
"fields":"temperature"
}
}
}
}
}
}``````
heat.json
``````// -*- mode: javascript -*-
{
"Name": "Thermo dynamics",
"ShortName":"ThermoDyn",
"Parameters":
{
"T_heater":310,
"T_ext":280
},
"Materials":
{
"air":
{
"markers":["air"],
"physics":["fluid","heat"],
"rho":"1",
"mu":"2.65e-2",
"k":"0.03",
"Cp":"1004",
"beta":"0.003660" //0.00006900
}
},
"BoundaryConditions":
{
"temperature":
{
"Dirichlet":
{
"heater1": { "expr":"T_heater:T_heater" }
},
"Robin":
{
"exterior-walls":
{
"expr1":"1.0/(0.06+0.01/0.5 + 0.3/0.8 + 0.20/0.032 +0.016/0.313 +0.14)",// h coeff
"expr2":"T_ext:T_ext"// temperature exterior
}
}
}
},
"PostProcess":
{
"Measures":
{
"Points":
{
"pointA":
{
"coord":"{0.05,0.25}",
"fields":"temperature"
},
"pointB":
{
"coord":"{0.55,0.25}",
"fields":"temperature"
},
"pointC":
{
"coord":"{1.05,0.25}",
"fields":"temperature"
},
"pointD":
{
"coord":"{1.55,0.25}",
"fields":"temperature"
},
"pointE":
{
"coord":"{0.05,0.5}",
"fields":"temperature"
},
"pointF":
{
"coord":"{0.55,0.5}",
"fields":"temperature"
},
"pointG":
{
"coord":"{1.05,0.5}",
"fields":"temperature"
},
"pointH":
{
"coord":"{1.55,0.5}",
"fields":"temperature"
},
"pointI":
{
"coord":"{0.05,0.75}",
"fields":"temperature"
},
"pointK":
{
"coord":"{0.55,0.75}",
"fields":"temperature"
},
"pointL":
{
"coord":"{1.05,0.75}",
"fields":"temperature"
},
"pointM":
{
"coord":"{1.55,0.75}",
"fields":"temperature"
}
}
}
}
}``````

See Case 1.

See Case 2.

See Case 3.

## 5. Bibliographie

• Validation expérimentale de modèles : application aux bâtiments basse consommation - Stephanie Bontemps

• Méthodes des bases réduites pour la qualité de l’air en milieu urbain - Janelle K.HAMMOND