Projet de 4fastsim (method of PBDW)

1. Introduction

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

\[\begin{equation} u^{bk}(\textbf{p}) + \eta \end{equation}\]

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

\[\begin{equation} u^{bk}_N(\textbf{p}) = \sum^{N}_{i=1}\beta_i(\textbf{p}) u^{bk}_N(\textbf{p}_i) \end{equation}\]

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->addParameterInModelProperties( param.first, param.second);

     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->addParameterInModelProperties( param.first, param.second );

         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 g = integrate(_range=elements( heat->mesh()), _expr= inner(gradv(u),gradv(v))).evaluate()(0,0);
   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" );
   pbdwdatafromaerooptions.add( toolboxes_options("heat-fluid") );
   pbdwdatafromaerooptions.add( toolboxes_options("heat") );
   pbdwdatafromaerooptions.add_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,
                    _about=about(_name="feelpp_ibat_pbdw_datafromaero",
                                 _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(gradt(qm),grad(v)));
     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) );
   double h1 = normH1(_range=elements(heat->mesh()), _expr=idv(u_true)-idv(u_MN), _grad_expr=gradv(u_true)-gradv(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";
   double c_h1 = normH1(_range=elements(heat->mesh()), _expr=idv(u_true), _grad_expr=gradv(u_true) );
   //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()),
       _expr=idv(u_true)-idv(ZN_orthn[i]) - idv(UM[j]), _grad_expr=gradv(u_true)-gradv(ZN_orthn[i]) - gradv(UM[j])));

       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 -> add("u_true", 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 -> add("u_err", u_err);
   er->save();


   auto d = exporter( _mesh=heat->mesh(),_name="u_h" );
   for(int i = 0; i < ZN.size(); i++)
       d->add( (boost::format("Z_%1%")%i).str(), ZN[i] );
   d->save();

   auto q = exporter( _mesh=heat->mesh(),_name="qm" );
   for(int i = 0; i < UM.size(); i++)
       q->add( (boost::format("q_%1%")%i).str(), UM[i] );
   q->save();

   auto n = exporter( _mesh=heat->mesh(),_name="Znorthn" );
   for(int i = 0; i < ZN_orthn.size(); i++)
       n->add( (boost::format("Znorthn_%1%")%i).str(), ZN_orthn[i] );
   n->save();

   auto e = exporter( _mesh=heat->mesh(), _name="u_MN");
   e -> add("solution PBDW", 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->add("qm_total", QM);
   f->save();


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

   auto u_velocity = heatFluid->fluidModel()->fieldVelocity();
   auto h = exporter( _mesh=heat->mesh(),_name="u_velocity" );
   h->add("u_velocity", 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