Accueil | Cast3M

* fichier :  ther9.dgibi
************************************************************************
************************************************************************
*************************************************************************
* NOM         : THER9
* DESCRIPTION : Cas-test du p-laplacien (p>=1) inspiré par \cite[chapitre
*               8]{saramito}.
*
*               Teste la thermique en massif 2D conduction anisotrope.
*
*
*               Il s'agit de l'équation non-linéaire suivante :
*
*               - div k(T) grad T = f
*
*               avec terme source : f = 1
*               et coefficient de conductivité : k(T) = |grad T|^(p-2)
*               condition aux limites T = 0 sur le bord du domaine.
*
*               Cette équation est intéressante à plusieurs titres :
*               - la régularité de la solution varie avec p ;
*               - pour p=2, on a le Laplacien normal
*               - pour p=1, on a l'opérateur courbure moyenne
*               - pour p-> +\inf, la solution T tend vers la fonction
*                 distance au bord \cite{belyaev}.
*               - suivant les valeurs de p, le comportement des
*                 algorithmes de résolution non-linéaire (point fixe,
*                 Newton, accélération de convergence) varie.
*
*               Pour la méthode de Newton, la matrice tangente du
*               p-Laplacien s'exprime sous la forme d'un Laplacien avec
*               un tenseur de diffusion anisotrope :
*
*               k2(T) = k(T) I + (p-2) |grad T|^(p-4) gradT \ptensoriel gradT
*
*
*
*               Sur la d-sphère unité, pour p>1, on a une solution
*               radiale analytique de l'équation :
*                             T = C (1-r^a)
*
*                     avec : a = p / (p-1)
*                            C = d^(-1/(p-1)) / a
*
*               En effet, pour une solution T radiale ne dépendant que
*               de r, le p-Laplacien s'exprime comme \cite{franzina} :
*               -\delta_p T(r) = -(p-1)   |T'(r)|^(p-2) T''(r)
*                                -(d-1)/r |T'(r)|^(p-2) T'(r)
*
*
*               Bibliographie (bibtex) :
*@unpublished{saramito,
*  TITLE = {{Efficient C++ finite element computing with Rheolef}},
*  AUTHOR = {Saramito, Pierre},
*  URL = {https://cel.archives-ouvertes.fr/cel-00573970},
*  NOTE = {Notes de cours de DEA accessibles depuis : \url{https://cel.archives-ouvertes.fr/cel-00573970} },
*  TYPE = {DEA},
*  ADDRESS = {Grenoble, France, France},
*  PAGES = {176},
*  YEAR = {2016},
*  MONTH = Jun,
*  KEYWORDS = {Finite Element Method FEM ; C++ LIBRARY ; Navier Stokes Equations ; Characteristic Method ; Elasticity equations},
*  PDF = {https://cel.archives-ouvertes.fr/cel-00573970/file/rheolef.pdf},
*  HAL_ID = {cel-00573970},
*  HAL_VERSION = {v13},
*}
*
*@article{franzina2010existence,
*  title={Existence and uniqueness for a p-laplacian nonlinear eigenvalue problem},
*  author={Franzina, Giovanni and Lamberti, Pier Domenico},
*  journal={Electron. J. Differential Equations},
*  volume={26},
*  number={10},
*  year={2010}
*}
*
*
*@article{belyaev,
*author = {Belyaev, Alexander G. and Fayolle, Pierre-Alain},
*title = {On Variational and PDE-Based Distance Function Approximations},
*journal = {Computer Graphics Forum},
*volume = {34},
*number = {8},
*issn = {1467-8659},
*url = {http://dx.doi.org/10.1111/cgf.12611},
*doi = {10.1111/cgf.12611},
*pages = {104--118},
*year = {2015},
*}
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SEMT/LTA)
*               mél : stephane.gounand@cea.fr
**********************************************************************
* VERSION    : v1, 02/05/2017, version initiale
* HISTORIQUE : v1, 02/05/2017, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
'OPTI' 'ECHO' 1 ;
'SAUT' 'PAGE' ;
interact = faux ;
graph    = faux ;
*
* Procédure d'accélération de la convergence.
*
'DEBP' ACCCVG ;
'ARGU' res*'CHPOINT' ;
'ARGU' tabacc*'TABLE' ;
'SI' ('NON' ('EXIS' tabacc 'iter')) ;
   tabacc . 'iter'   = 0 ;
   tabacc . 'znacce' = 4 ;
   tabacc . 'itdep'  = 3 ;
   tabacc . 'acfp1'  = 'COPI' ('*' res 0.) ;
   tabacc . 'acfp2'  = tabacc . 'acfp1' ;
   tabacc . 'acfp3'  = tabacc . 'acfp1' ;
   tabacc . 'acfep1'  = tabacc . 'acfp1' ;
   tabacc . 'acfep2'  = tabacc . 'acfp1' ;
   tabacc . 'correc'  = 0. ;
   tabacc . 'freap'   = 0. ;
'FINS' ;   
it = '+' (tabacc . 'iter') 1 ;
znacce = tabacc . 'znacce' ;
itdep  = tabacc . 'itdep'  ;
tabacc . 'iter' = it ;
tabacc . 'correcp'  = tabacc . 'correc' ;
tabacc . 'correc'   = 0. ;
tabacc . 'acfp0'  = 'ENLEVER' ('-' res (tabacc . 'freap')) 'FLX' ;
tabacc . 'acfep0' = '-' (tabacc . 'acfp0') (tabacc . 'correcp') ;
'SI' ('ET' ('MULT' it znacce) ('>' it itdep)) ;
   'MESS' 'acceleration' ;
   correc = 'ACT3' (tabacc . 'acfep2')  (tabacc . 'acfep1')
                   (tabacc . 'acfep0')  (tabacc . 'acfp3')                    
                   (tabacc . 'acfp2')   (tabacc . 'acfp1')
                   (tabacc . 'acfp0')   ;
   res = '-' res correc ;
'FINS' ;
'SI' ('>' it 3) ;
   'DETR' (tabacc . 'acfp3') ;
   'DETR' (tabacc . 'acfep2') ;
'FINSI' ;   
tabacc . 'acfp3'  = tabacc . 'acfp2 ' ;
tabacc . 'acfp2'  = tabacc . 'acfp1 ' ;
tabacc . 'acfp1'  = tabacc . 'acfp0 ' ;
tabacc . 'acfep2'  = tabacc . 'acfep1 ' ;
tabacc . 'acfep1'  = tabacc . 'acfep0 ' ;
'RESP' res tabacc ;
'FINP' ;
*
* Options globales
*
'OPTI' 'DIME' 2 ;
'OPTI' 'ELEM' 'TRI3' ;
*
* Maillage
*
* nx : nombre de maille par côté
* idom : type de domaine 0 (carré) ou 1 (quart de cercle)
*        Sur 1, il y a une solution analytique pour p>1
nx = 40 ;
idom = 1 ;
* Physique
* vp : exposant p du p-laplacien
vp = 3. ;
* Numérique
* nitmax : nombre d'itérations max.
* tol   : tolérance pour la convergence de la boucle non-linéaire
* omega : facteur de relaxation
* iktan : type de matrice tangente :
*         0 laplacien ; 1 point fixe ou ; 2 newton
nitmax = 100 ;
tol = 1.D-10 ;
*omega = '/' 2. vp ;
*iktan = 1 ;
omega = 1. ;
iktan = 2 ;
lacc  = 0 ; tabacc = 'TABL' ;
 
chpar = 'CHAI' 'p=' vp ' omega=' omega ' iktan=' iktan ' lacc=' lacc ; 
 
'MESS' '*************************' ;
'MESS' '* P laplacien ' chpar ;
'MESS' '*************************' ;
* Points
'OPTI' 'DENS' ('/' 1. ('FLOT' nx)) ;
pA = 0. 0. ;
pB = 1. 0. ;
pC = 1. 1. ;
pD = 0. 1. ;
*
* Maillage carré 
*
'SI' ('EGA' idom 0) ;
* Lignes
   lAB = 'DROI' nx pA pB ;
   lBC = 'DROI' nx pB pC ;
   lCD = 'DROI' nx pC pD ;
   lDA = 'DROI' nx pD pA ;
* Surface
*mt = 'DALL' lAB lBC lCD lDA ;
   mt = 'SURF' (lAB 'ET' lBC 'ET' lCD 'ET' lDA) ;
* Blocages   
   mailblo = 'CONT' mt ;
'FINS' ;
*
* Maillage quart de cercle
*
'SI' ('EGA' idom 1) ;
* Lignes
   lAB = 'DROI' pA pB ;
   lBD = 'CERC' pB pA pD ;
   lDA = 'DROI' pD pA ;
* Surface
*mt = 'DALL' lAB lBC lCD lDA ;
   mt = 'SURF' (lAB 'ET' lBD 'ET' lDA) ;
   mailblo = lBD ;
'FINS' ;   
*
'SI' graph ;
   'TRAC' mt 'TITR' 'MAILLAGE' ;
'FINS' ;   
*
* Modele et materiau initial
*
modt  = 'MODE' mt 'THERMIQUE' 'ISOTROPE' ;
modt2 = 'MODE' mt 'THERMIQUE' 'ANISOTROPE' 'CONS' 'T2';
matt = 'MATE' modt 'K' 1. ;
*
mlap0 = 'COND' modt matt ;
* terme source
tsou = 'SOUR' modt 1. mt ;
* conditions aux limites
mcli = 'BLOQ' 'T' mailblo ;
* résolution
mat = mlap0 'ET' mcli ;
smb = tsou ;
sol = 'RESO' mat smb ;
solt = 'EXCO' 'T' sol 'T' ;
*ch0 =  'Solution au pas 0' ;
*'TRAC' solt mt 'TITR' ch0 ;
ipas = 0 ; soli = sol ; 
lnres = 'PROG' ;
cgtx = 'MOTS' 'T,X' ; cgty = 'MOTS' 'T,Y' ;
cgt  = cgtx 'ET' cgty ;
ck11 = 'MOTS' 'K11' ; ck21 = 'MOTS' 'K21' ; ck22 = 'MOTS' 'K22' ;
*
* Boucle pour la non-linearite
*
'REPE' bcl nitmax ;
   gsoli = 'GRAD' modt soli ;
   ngsoli = 'PSCA' gsoli gsoli cgt cgt ;
* Nouveau matériau avec la nouvelle conductivité
   chk = ngsoli '**' ('/' ('-' vp 2.) 2.) ;
   matt = 'MATE' modt 'K' chk ;
* Matrice et calcul du résidu 
   mlap = 'COND' modt matt ;
   mat =  mlap 'ET' mcli ;
   res = tsou '-' ('*' mat soli) ;
*   
   nres = 'MAXI' res 'ABS' 'AVEC' ('MOTS' 'Q') ;
   lnres = 'ET' lnres ('PROG' nres) ;
   'SI' ('<' nres tol) ;
      'MESS' ('CHAI' 'Pas i=' ipas ' max|Ri|=' nres ' < tol=' tol) ;
      'QUIT' bcl ;
   'SINO' ;      
      'MESS' ('CHAI' 'Pas i=' ipas ' max|Ri|=' nres) ;
   'FINS' ;
* Matrice tangente = rigidite au pas 0
   'SI' ('EGA' iktan 0) ;
      mattan = mlap0 ;
   'FINS' ;
* Matrice tangente = Matrice de point fixe
   'SI' ('EGA' iktan 1) ;
      mattan = mlap ;
   'FINS' ;
   'SI' ('EGA' iktan 2) ;
* Matrice tangente = Matrice de Newton
* 2 \eta' (|grad T|^2) (grad T \ptens grad T)
*  avec \eta'(z) = \frac{p-2}{2} z^{\frac{p-4}{2}}
      coeff = '*' (ngsoli '**' ('/' ('-' vp 4.) 2.))
                  ('-' vp 2.) ;
* On change le type de gsoli sinon il veut une syntaxe avec un modele
      gsoli = 'CHAN' 'TYPE' gsoli ' ' ;
      vk11  = '*' gsoli gsoli cgtx cgtx ck11 ;
      vk21  = '*' gsoli gsoli cgty cgtx ck21 ;
      vk22  = '*' gsoli gsoli cgty cgty ck22 ;
      matt2 = 'MATE' modt2 'DIRECTION' (1. 0.) 'PARALLELE'
                     'K11' ('CHAN' ('*' vk11 coeff) 'CONS' 'T2')
                     'K21' ('CHAN' ('*' vk21 coeff) 'CONS' 'T2')
                     'K22' ('CHAN' ('*' vk22 coeff) 'CONS' 'T2') ;
      matnew2 = 'COND' modt2 matt2 ;
      mattan = mlap 'ET' matnew2 ;
   'FINS' ;      
   mattot = ('*' mattan ('/' 1. omega)) 'ET' mcli ;
   smb = res ;
   'SI' ('EGA' lacc 1) ; smb tabacc = ACCCVG smb tabacc ; 'FINS' ;
   dsol = 'RESO' mattot smb ;
   'SI' ('EGA' lacc 1) ; tabacc . 'freap' = '*' mcli dsol ; 'FINS' ;
   solip = '+' soli dsol ;
*   solipt = 'EXCO' 'T' solip 'T' ;
*   chip =  'CHAI' 'Solution  vp=' vp ' au pas ' ('+' ipas 1) ;
*   'TRAC' solipt mt 'TITR' chip 'NCLK' ;
   soli = solip ;
   ipas = '+' ipas 1 ;
'FIN' bcl ;
*
* Courbe du résidu
*
npas = 'DIME' lnres ;
'SI' ('>' npas 1) ;
   lpas = 'PROG' 0. 'PAS' 1. 'NPAS' ('-' npas 1) ;
   llnres = '/' ('LOG' lnres) ('LOG' 10.) ;
   evl = 'EVOL' 'MANU' lpas llnres ;
   'SI' graph ;
      titx = 'iter' ; tity = 'Log. |res|' ;
      tit = 'CHAI' tity '(' titx ')' ' ' chpar ;
      'DESS' evl 'TITR' tit 'TITX' titx 'TITY' tity ;
   'FINS' ;      
'FINS' ;   
*
* Tracé en élévation
*
'SI' graph ;
   'OPTI' 'DIME' 3 ;
   solit = 'EXCO' 'T' soli ;                   
   depz = 'NOMC' solit 'UZ' ;
   'FORM' ('/' depz ('MAXI' depz 'ABS')) ;
   'OPTI' 'OEIL' (20. 30. 10.) ;
   'OPTI' 'ISOV' 'SULI' ;
   chi =  'CHAI' 'Solution  vp=' vp ' au pas ' ipas ;
   'TRAC' solit mt ('CONT' mt) 30 'TITR' chi ;
   'OPTI' 'DIME' 2 ;
'FINS' ;
*
* Comparaison à la solution analytique
*
'SI' ('EGA' idom 1) ;
   x y = 'COOR' mt ;
   r2 = (x '*' x) '+' (y '*' y) ;
   expo  = vp '/' ('-' vp 1.) ;
   vd = 'VALE' 'DIME' ;
   Const = '/' (vd '**' ('/' -1 ('-' vp 1.)))
               expo ;
   solana = '*' ('-' 1. ('**' r2 ('/' expo 2.)))
                const ;
   solnum = 'EXCO' 'T' soli ;                   
   dsol = '-' solana solnum ;
   'SI' graph ;
      'TRAC' dsol mt ('CONT' mt) 30
         'TITR' 'Difference analytique-numerique' ;
   'FINS' ;         
   mdsol = 'MAXI' dsol 'ABS' ;
   ch = 'CHAI' 'Maxi. |analytique - numerique|=' mdsol ' ' chpar;
   'MESS' ch ;
'FINS' ;
*
* Tests :
* 1) La méthode de Newton converge en 7-8 itérations pour p=3 (cf. Saramito
* p.137)  
* 2) L'erreur par rapport à la solution analytique vaut environ 4.e-4
*    avec 40 mailles
*
lok = vrai ;
lok = 'ET' lok (npas '&lt;EG' 8) ;
'SI' ('NON' lok) ;
   ch = 'CHAI' 'Newton : unexpectedly slow convergence' ;
   'ERRE' ch ;
'FINS' ;
lok = 'ET' lok (mdsol '<' 5.e-4) ;
'SI' ('NON' lok) ;
   ch = 'CHAI' 'Comparison with analytical solution : bad results' ;
   'ERRE' ch ;
'FINS' ;
'SI' ('NON' lok) ;
   'ERREUR' 5 ;
'SINON' ;   
   'MESSAGE' ('CHAINE' 'Test ok !') ;
'FINSI' ;
*
'SI' interact ;
   'OPTION' 'DONN' 5 'ECHO' 1 ;
'FINSI' ;
*
* End of dgibi file THER9
*
'FIN' ;