test/test6905/test914_Xternal_4_Variables.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2020 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 #===============================================================================
  24 import numpy
  25 # donnees numeriques du modele
  26 xdeb = 0.
  27 xfin = 1.
  28
  29 nx = 500
  30 dx=(xfin-xdeb)/nx
  31 dim = nx+1
  32 p = 10 # nb de mesures
  33
  34 class EqChaleur(object):
  35
  36     def __init__(self,alpha):
  37
  38         self.alpha = float(alpha)
  39
  40         idx2 = 1./dx**2
  41         self.A=2*idx2*numpy.identity(dim)
  42         for i in range(1,dim-2):
  43             self.A[i,i+1]=-idx2
  44             self.A[i+1,i]=-idx2
  45
  46     def solve(self,S):
  47         return numpy.linalg.solve(self.alpha*self.A,S)
  48
  49     def get2ndMbre(self):
  50
  51         S = numpy.zeros((dim,))
  52         c = numpy.pi*dx
  53         pi2 = numpy.pi**2
  54
  55         for i in range(dim):
  56             S[i] = pi2*numpy.sin(i*c)
  57         return S
  58
  59     def extract(self,champs,nbobs):
  60
  61         incr = int(nx/nbobs)
  62         return champs[::incr]
  63
  64     def obs(self,nbobs):
  65
  66         X = numpy.linspace(xdeb,xfin,dim)
  67         S = self.get2ndMbre()
  68         T = self.solve(S)
  69         return self.extract(T,nbobs)
  70
  71     def disp(self):
  72
  73         X = numpy.linspace(xdeb,xfin,dim)
  74         S = self.get2ndMbre()
  75         T = self.solve(S)
  76
  77         from matplotlib import pylab
  78         pylab.plot(X,T)
  79         pylab.show()
  80
  81 class ExempleOperateurChaleur:
  82     """
  83     Modelisation d'operateur non lineaire pour la diffusion de la chaleur dans
  84     un mur d'epaisseur e, avec une production de chaleur r, pour un
  85     environnement a temperature Ts. Les equations sont :
  86
  87         T = Ts - r/2k x^2 + r e /2k  x
  88     """
  89     def __init__(self, e = 10., Ts = 100., nbpt = 1000, nbobs = 10):
  90         self.e     = float(e)
  91         self.Ts    = float(Ts)
  92         self.nbpt  = int(nbpt)
  93         self.dx    = self.e/(self.nbpt-1)
  94         self.nbobs = min(int(nbobs),self.nbpt)
  95
  96     def Direct(self, XX ):
  97         if type(XX) is type(numpy.matrix([])):  alpha = XX.A1
  98         elif type(XX) is type(numpy.array([])): alpha = numpy.matrix(XX).A1
  99         else:                                alpha = XX
 100         #
 101         eq=EqChaleur(alpha)
 102         HX = eq.obs(self.nbobs)
 103         #
 104         return numpy.array( HX )
 105
 106     def DirectExp(self, XX ): # Version en Log/Exp
 107         if type(XX) is type(numpy.matrix([])):  alpha = XX.A1
 108         elif type(XX) is type(numpy.array([])): alpha = numpy.matrix(XX).A1
 109         else:                                   alpha = XX
 110         alpha = numpy.exp(alpha)
 111         #
 112         eq=EqChaleur(alpha)
 113         HX = eq.obs(self.nbobs)
 114         #
 115         return numpy.array( HX )
 116
 117
 118     def Analytique(self, XX ):
 119         amort = float(XX)
 120         listresult = []
 121         for i in range(self.nbobs+1) :
 122             ex = i * self.e/(self.nbobs)
 123             eT = self.Ts - amort * ex**2 + self.e*amort*ex
 124             listresult.append(eT)
 125         return numpy.array(listresult)
 126
 127     def Tangent(self, paire ):
 128         (X, dX) = paire
 129         pass
 130
 131     def Adjoint(self, paire ):
 132         (X, Y) = paire
 133         pass
 134
 135 #===============================================================================
 136
 137 CheckingPoint = [2.]
 138
 139 AlgorithmParameters = {
 140     "EpsilonMinimumExponent" : -10,
 141     # "PlotAndSave" : True,
 142     "SetSeed" : 1000,
 143     "InitialDirection":CheckingPoint,
 144     "AmplitudeOfInitialDirection":0.5,
 145     }
 146
 147 taille = 1000
 148 OP = ExempleOperateurChaleur(nbpt = taille)
 149
 150 DirectOperator = OP.Direct
 151
 152 # JPA : probleme : il faut que le nom d'origine, qui est une méthode de classe,
 153 # existe explicitement pour permettre le parallelisme !
 154 Direct = DirectOperator