Cases/crue-ADAO-3DVAR-1parametre.py

   1 # -*- coding: utf-8 -*-\r
   2 # Copyright (C) 2019 - Jean-Philippe Argaud\r
   3 \r
   4 import numpy as np\r
   5 import unittest\r
   6 np.set_printoptions(precision=3)\r
   7 from adao import adaoBuilder\r
   8 \r
   9 # Les entrées\r
  10 nbobs = 4\r
  11 Qobs = np.array([10.,20.,30.,40.])\r
  12 Hobs = np.array([0.19694513, 0.298513, 0.38073079, 0.45246109])\r
  13 \r
  14 # Définit les paramètres de l'AD\r
  15 KsInitial = 20.\r
  16 ZvInitial = 49.\r
  17 ZmInitial = 51.\r
  18 thetaB = [KsInitial,]\r
  19 \r
  20 def functionCrue(Q, K_s):\r
  21     L = 5.0e3\r
  22     B = 300.0\r
  23     Z_v = ZvInitial\r
  24     Z_m = ZmInitial\r
  25     alpha = (Z_m - Z_v)/L\r
  26     H = (Q/(K_s*B*np.sqrt(alpha)))**(3.0/5.0)\r
  27     return H\r
  28 \r
  29 def obsFunction(theta):\r
  30     # Evaluation de la sortie du modèle\r
  31     K_s = float(theta)\r
  32     sortie = np.zeros((nbobs,))\r
  33     for i in range(nbobs):\r
  34         sortie[i] = functionCrue(Qobs[i],K_s)\r
  35     return sortie\r
  36 \r
  37 def obsFunctionMulti(thetas):\r
  38     ret = []\r
  39     for elt in thetas:\r
  40         ret.append( np.array(obsFunction(elt)) )\r
  41     return ret\r
  42 \r
  43 class TestRefFloodTest(unittest.TestCase):\r
  44 \r
  45     def test0(self):\r
  46         # define the problem bounds\r
  47         boundsMin = [20.,]\r
  48         boundsMax = [40.,]\r
  49 \r
  50         # Ecart-type des paramètres\r
  51         sigmaKs = 5.e10\r
  52         sigmaZv = 1.\r
  53         sigmaZm = 1.\r
  54         sigmaTheta = [sigmaKs,]\r
  55 \r
  56         # Ecart-type des observations\r
  57         sigmaH = 0.5 # (m^2)\r
  58 \r
  59         case = adaoBuilder.New()\r
  60         case.set( 'AlgorithmParameters', Algorithm='3DVAR' )\r
  61         case.set( 'Background',          Vector=thetaB)\r
  62         case.set( 'BackgroundError',     DiagonalSparseMatrix=sigmaTheta )\r
  63         case.set( 'Observation',         Vector=Hobs)\r
  64         case.set( 'ObservationError',    ScalarSparseMatrix=sigmaH )\r
  65         case.set( 'ObservationOperator', OneFunction= obsFunctionMulti, InputFunctionAsMulti = True)\r
  66         case.execute()\r
  67         res = case.get("Analysis")[-1][0]\r
  68         self.assertAlmostEqual(res,25.,6)\r
  69         pass\r
  70     pass\r
  71 \r
  72 if __name__ == "__main__":\r
  73     unittest.main()\r
  74     pass\r
  75 \r