doc/doxygen/doxfiles/interpolation/interptheory.dox

   1 /*!
   2 \page InterpKerRemapGlobal Linear remapping
   3
   4 For fields with polynomial representation on each cell, the components of the discretized field  \f$ \phi_s \f$ on the source side can be expressed as linear combinations of the components of the discretized field \f$ \phi_t \f$ on the target side, in terms of a matrix-vector product:
   5
   6 \f[
   7  \phi_t=W.\phi_s.
   8 \f]
   9
  10 \f$W\f$ is called the \anchor interpolationmatrix interpolation matrix.
  11 The objective of interpolators is to compute the matrix W depending on their physical
  12 properties (\ref IntExtFields) and their mesh discretization (on cells P0, on nodes P1,...).
  13
  14 \section ConsInterp Conservative interpolation
  15
  16 At the basis of many CFD numerical schemes is the fact that physical
  17 quantities such as density, momentum per unit volume or energy per
  18 unit volume obey some balance laws that should be preserved at the
  19 discrete level on every cell.
  20
  21 It is therefore often desired that the process interpolation preserve the
  22 integral of \f$ \phi \f$ on any domain. At the discrete level, for any
  23 target cell \f$ T_i \f$, the following \b general \b interpolation \b
  24 formula has to be satisfied :
  25 \anchor InterpKerGenralEq
  26 \f[
  27 \int_{T_i} \phi_t = \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi_s.
  28 \f]
  29
  30 This equation is used to compute \f$ W_{ij} \f$, based on the fields representation ( P0, P1, P1d etc..) and the
  31 geometry of source and target mesh cells.
  32
  33 \section MeshOverlap Mesh overlapping
  34
  35 Another important property of the interpolation process is the maximum principle: the field values resulting from the interpolation should remain between the upper and lower bounds of the original field.
  36 When interpolation is performed between a source mesh S and a target
  37 mesh T the aspect of overlapping is important. In fact if any cell of
  38 of S is fully overlapped by cells of T and inversely any cell of T is
  39 fully overlapped by cells of S that is
  40 \f[
  41 \sum_{S_j} Vol(T_i\cap S_j) = Vol(T_i),\hspace{1cm} and \hspace{1cm} \sum_{T_i} Vol(S_j\cap T_i) = Vol(S_j)
  42 \f]
  43 then the meshes S and T are said to be \b
  44 overlapping. In this case the two formulas in a given column in the table below give the same
  45 result. All intensive formulas result in the same output, and all the extensive formulas give also the same output.
  46
  47 The ideal interpolation algorithm should be conservative and respect the maximum principle. However such an algorithm can be impossible to design if the two meshes do not overlap. When the meshes do not overlap, using either \f$Vol(T_i)\f$ or \f$\sum_{S_j} Vol(T_i\cap S_j)\f$ one obtains an algorithm that respects either the conservativity or the maximum principle (see the nature of field \ref TableNatureOfField "summary table").
  48
  49
  50 \section InterpKerRemapInt Linear conservative remapping of P0 (cell based) fields
  51
  52 We assume that the field is represented by a vector with a discrete value on each cell.
  53 This value can represent either
  54 - an average value of the field in the cell (average density, velocity or temperature in the cell) in which case the representation is said to be \b intensive,
  55 - an integrated value over the cell (total mass, power of the cell) in which case the representation is said to be \b extensive
  56
  57 \section InterpKerP0P0Int cell-cell (P0->P0) conservative remapping of intensive fields
  58
  59 For intensive fields such as mass density or power density, the
  60 left hand side in the \ref InterpKerGenralEq "general interpolation equation" becomes :
  61
  62 \f[
  63 \int_{T_i} \phi = Vol(T_i).\phi_{T_i}.
  64 \f]
  65
  66 Here Vol represents the volume when the mesh dimension is equal to 3, the
  67 area when mesh dimension is equal to 2, and length when mesh dimension is equal to 1.
  68
  69 In the \ref InterpKerGenralEq "general interpolation equation" the
  70 right hand side becomes :
  71
  72 \f[
  73 \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi = \sum_{S_j\cap T_i \neq \emptyset} {Vol(T_i\cap S_j)}.\phi_{S_j}.
  74 \f]
  75
  76 As the field values are constant on each
  77 cell, the coefficients of the linear remapping matrix \f$ W \f$ are
  78 given by the formula :
  79
  80 \f[
  81  W_{ij}=\frac{Vol(T_i\cap S_j)}{ Vol(T_i) }.
  82 \f]
  83
  84
  85 \section InterpKerP0P0Ext cell-cell (P0->P0) conservative remapping of extensive fields
  86
  87 In code coupling from neutronics to hydraulics, \b extensive field
  88 of power is exchanged and the total power should remain the same.
  89 The discrete values of the field represent the total power contained in the cell.
  90 Hence in the \ref InterpKerGenralEq "general interpolation equation" the
  91 left hand side becomes :
  92
  93 \f[
  94 \int_{T_i} \phi = P_{T_i},
  95 \f]
  96
  97 while the right hand side is now :
  98
  99 \f[
 100 \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi =
 101 \sum_{S_j\cap T_i \neq \emptyset} \frac{Vol(T_i\cap S_j)}{ Vol(S_j)}.P_{S_j}.
 102 \f]
 103
 104 The coefficients of the linear remapping matrix \f$ W \f$ are then
 105 given by the formula :
 106
 107 \f[
 108  W_{ij}=\frac{Vol(T_i\cap S_j)}{  Vol(S_j) }.
 109 \f]
 110
 111 \section TableNatureOfField Summary
 112 In the case of fields with a P0 representation (cell based) and when the meshes do not overlap, the scheme is either conservative or maximum preserving (not both). Depending on the \ref NatureOfField the interpolation coefficients take the following value:
 113
 114  * <TABLE BORDER=1 >
 115  * <TR><TD> </TD><TD>Intensive</TD><TD> extensive </TD></TR>
 116  * <TR><TD> Conservation</TD><TD> \f[\frac{Vol(T_i\cap S_j)}{ Vol(T_i)}\f] <br /> \ref TableNatureOfFieldExampleRevIntegral "RevIntegral" </TD><TD> \f[ \frac{Vol(T_i\cap S_j)}{ \sum_{T_i} Vol(S_j\cap T_i) }\f] <br /> \ref TableNatureOfFieldExampleIntegralGlobConstraint "IntegralGlobConstraint" </TD></TR>
 117  * <TR><TD> Maximum principle </TD><TD> \f[\frac{Vol(T_i\cap S_j)}{ \sum_{S_j} Vol(T_i\cap S_j)}\f] <br /> \ref TableNatureOfFieldExampleConservVol "ConservativeVolumic" </TD><TD>  \f[\frac{Vol(T_i\cap S_j)}{  Vol(S_j) }\f] <br /> \ref TableNatureOfFieldExampleIntegral "Integral"</TD></TR>
 118  *</TABLE>
 119
 120 \section TableNatureOfFieldExample Illustration of a non overlapping P0P0 interpolation
 121
 122 Let's consider the following case with a source mesh containing two cells and a target mesh containing one cell.
 123 Let's consider a field FS on cells on the source mesh that we want to interpolate on the target mesh.
 124
 125 The value of FS on the cell#0 is 4 and the value on the cell#1 is 100.
 126
 127 The aim here is to compute the interpolated field FT on the target mesh of field FS depending on the \ref NatureOfField "nature of the field".
 128
 129 \anchor TableNatureOfFieldEx1
 130 \image html NonOverlapping.png "An example of non overlapping intersection of two meshes."
 131
 132 The first step of the interpolation leads to the following M1 matrix :
 133
 134 \f[
 135     M1=\left[\begin{tabular}{cc}
 136     0.125 & 0.75 \\
 137     \end{tabular}\right]
 138     \f]
 139
 140 \subsection TableNatureOfFieldExampleConservVol Conservative volumic case
 141
 142 If we apply the formula \ref TableNatureOfField "above" it leads to the following \f$ M_{Conservative Volumic} \f$ matrix :
 143
 144 \f[
 145     M_{Conservative Volumic}=\left[\begin{tabular}{cc}
 146     $\displaystyle{\frac{0.125}{0.125+0.75}}$ &
 147     $\displaystyle{\frac{0.75}{0.125+0.75}}$ \\
 148     \end{tabular}\right]=\left[\begin{tabular}{cc}
 149     0.14286 & 0.85714 \\
 150     \end{tabular}\right]
 151 \f]
 152 \f[
 153     FT=\left[\begin{tabular}{cc}
 154     $\displaystyle\frac{0.125}{0.875}$ & $\displaystyle\frac{0.75}{0.875}$ \\
 155     \end{tabular}\right].\left[\begin{tabular}{c}
 156     4 \\
 157     100 \\
 158     \end{tabular}\right]
 159     =\left[\begin{tabular}{c}
 160     86.28571\\
 161     \end{tabular}\right]
 162 \f]
 163
 164 As we can see here the maximum principle is respected.This nature of field is particularly recommended to interpolate an intensive
 165 field such as \b temperature or \b pressure.
 166
 167 \subsection TableNatureOfFieldExampleIntegral Integral case
 168
 169 If we apply the formula \ref TableNatureOfField "above" it leads to the following \f$ M_{Integral} \f$ matrix :
 170
 171 \f[
 172     M_{Integral}=\left[\begin{tabular}{cc}
 173     $\displaystyle{\frac{0.125}{9}}$ & $\displaystyle{\frac{0.75}{3}}$ \\
 174     \end{tabular}\right]=\left[\begin{tabular}{cc}
 175     0.013888 & 0.25 \\
 176     \end{tabular}\right]
 177 \f]
 178 \f[
 179     FT=\left[\begin{tabular}{cc}
 180     $\displaystyle{\frac{0.125}{9}}$ & $\displaystyle{\frac{0.75}{3}}$ \\
 181     \end{tabular}\right].\left[\begin{tabular}{c}
 182     4 \\
 183     100 \\
 184     \end{tabular}\right]
 185     =\left[\begin{tabular}{c}
 186     25.055\\
 187     \end{tabular}\right]
 188 \f]
 189
 190 This type of interpolation is typically recommended for the interpolation of \b power (\b NOT \b power \b density !) for
 191 a user who wants to conserve the quantity \b only on the intersecting part of the source mesh (the green part on the \ref TableNatureOfFieldEx1 "example")
 192
 193 This type of interpolation is equivalent to the computation of \f$ FS_{vol} \f$ followed by a multiplication by \f$ M1 \f$ where \f$ FS_{vol} \f$ is given by :
 194
 195 \f[
 196    FS_{vol}=\left[\begin{tabular}{c}
 197     $\displaystyle{\frac{4}{9}}$ \\
 198     $\displaystyle{\frac{100}{3}}$ \\
 199     \end{tabular}\right]
 200 \f]
 201
 202 In the particular case treated \ref TableNatureOfFieldEx1 "here", it means that only a power of 25.055 W is intercepted by the target cell !
 203
 204 So from the 104 W of the source field \f$ FS \f$, only 25.055 W are transmitted in the target field using this nature of field.
 205 In order to treat differently a power field, another policy, \ref TableNatureOfFieldExampleIntegralGlobConstraint "integral global constraint nature" is available.
 206
 207 \subsection TableNatureOfFieldExampleIntegralGlobConstraint Integral with global constraints case
 208
 209 If we apply the formula \ref TableNatureOfField "above" it leads to the following \f$ M_{IntegralGlobConstraint} \f$ matrix :
 210
 211 \f[
 212     M_{IntegralGlobConstraint}=\left[\begin{tabular}{cc}
 213     $\displaystyle{\frac{0.125}{0.125}}$ & ${\displaystyle\frac{0.75}{0.75}}$ \\
 214     \end{tabular}\right]=\left[\begin{tabular}{cc}
 215     1 & 1 \\
 216     \end{tabular}\right]
 217 \f]
 218 \f[
 219     FT=\left[\begin{tabular}{cc}
 220     1 & 1 \\
 221     \end{tabular}\right].\left[\begin{tabular}{c}
 222     4 \\
 223     100 \\
 224     \end{tabular}\right]
 225     =\left[\begin{tabular}{c}
 226     104\\
 227     \end{tabular}\right]
 228 \f]
 229
 230 This type of interpolation is typically recommended for the interpolation of \b power (\b NOT \b power \b density !) for
 231 a user who wants to \b conserve \b all \b the \b power in its source field. Here we have 104 W in source field, we have 104 W too,
 232 in the output target interpolated field.
 233
 234 \b BUT, As we can see here, the maximum principle is \b not respected here, because the target cell #0 has a value higher than the two
 235 intercepted source cells.
 236
 237 \subsection TableNatureOfFieldExampleRevIntegral Reverse integral case
 238
 239 If we apply the formula \ref TableNatureOfField "above" it leads to the following \f$ M_{RevIntegral} \f$ matrix :
 240
 241 \f[
 242     M_{RevIntegral}=\left[\begin{tabular}{cc}
 243     $\displaystyle{\frac{0.125}{1.5}}$ & $\displaystyle{\frac{0.75}{1.5}}$ \\
 244     \end{tabular}\right]=\left[\begin{tabular}{cc}
 245     0.083333 & 0.5 \\
 246     \end{tabular}\right]
 247 \f]
 248 \f[
 249     FT=\left[\begin{tabular}{cc}
 250     $\displaystyle{\frac{0.125}{1.5}}$ & $\displaystyle{\frac{0.75}{1.5}}$ \\
 251     \end{tabular}\right].\left[\begin{tabular}{c}
 252     4 \\
 253     100 \\
 254     \end{tabular}\right]
 255     =\left[\begin{tabular}{c}
 256     50.333\\
 257     \end{tabular}\right]
 258 \f]
 259
 260 This type of nature is particulary recommended to interpolate an intensive \b density
 261 field (moderator density, power density).
 262 The difference with \ref TableNatureOfFieldExampleConservVol "conservative volumic" seen above is that here the
 263 target field is homogenized to the \b whole target cell. It explains why this nature of field does not follow the maximum principle.
 264
 265 To illustrate the case, let's consider that \f$ FS \f$ is a power density field in \f$ W/m^2 \f$.
 266 With this nature of field the target cell #0 accumulates 0.125*4=0.5 W of power from the source cell #0 and 0.75*100=75 W of power from the source cell #1.
 267 It leads to 75.5 W of power on the \b whole target cell #0. So, the final power density is equal to 75.5/1.5=50.333 W/m^2.
 268
 269 */
 270
 271