Télécharger ther9.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : ther9.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *************************************************************************
  5. * NOM : THER9
  6. * DESCRIPTION : Cas-test du p-laplacien (p>=1) inspiré par \cite[chapitre
  7. * 8]{saramito}.
  8. *
  9. * Teste la thermique en massif 2D conduction anisotrope.
  10. *
  11. *
  12. * Il s'agit de l'équation non-linéaire suivante :
  13. *
  14. * - div k(T) grad T = f
  15. *
  16. * avec terme source : f = 1
  17. * et coefficient de conductivité : k(T) = |grad T|^(p-2)
  18. * condition aux limites T = 0 sur le bord du domaine.
  19. *
  20. * Cette équation est intéressante à plusieurs titres :
  21. * - la régularité de la solution varie avec p ;
  22. * - pour p=2, on a le Laplacien normal
  23. * - pour p=1, on a l'opérateur courbure moyenne
  24. * - pour p-> +\inf, la solution T tend vers la fonction
  25. * distance au bord \cite{belyaev}.
  26. * - suivant les valeurs de p, le comportement des
  27. * algorithmes de résolution non-linéaire (point fixe,
  28. * Newton, accélération de convergence) varie.
  29. *
  30. * Pour la méthode de Newton, la matrice tangente du
  31. * p-Laplacien s'exprime sous la forme d'un Laplacien avec
  32. * un tenseur de diffusion anisotrope :
  33. *
  34. * k2(T) = k(T) I + (p-2) |grad T|^(p-4) gradT \ptensoriel gradT
  35. *
  36. *
  37. *
  38. * Sur la d-sphère unité, pour p>1, on a une solution
  39. * radiale analytique de l'équation :
  40. * T = C (1-r^a)
  41. *
  42. * avec : a = p / (p-1)
  43. * C = d^(-1/(p-1)) / a
  44. *
  45. * En effet, pour une solution T radiale ne dépendant que
  46. * de r, le p-Laplacien s'exprime comme \cite{franzina} :
  47. * -\delta_p T(r) = -(p-1) |T'(r)|^(p-2) T''(r)
  48. * -(d-1)/r |T'(r)|^(p-2) T'(r)
  49. *
  50. *
  51. * Bibliographie (bibtex) :
  52. *@unpublished{saramito,
  53. * TITLE = {{Efficient C++ finite element computing with Rheolef}},
  54. * AUTHOR = {Saramito, Pierre},
  55. * URL = {https://cel.archives-ouvertes.fr/cel-00573970},
  56. * NOTE = {Notes de cours de DEA accessibles depuis : \url{https://cel.archives-ouvertes.fr/cel-00573970} },
  57. * TYPE = {DEA},
  58. * ADDRESS = {Grenoble, France, France},
  59. * PAGES = {176},
  60. * YEAR = {2016},
  61. * MONTH = Jun,
  62. * KEYWORDS = {Finite Element Method FEM ; C++ LIBRARY ; Navier Stokes Equations ; Characteristic Method ; Elasticity equations},
  63. * PDF = {https://cel.archives-ouvertes.fr/cel-00573970/file/rheolef.pdf},
  64. * HAL_ID = {cel-00573970},
  65. * HAL_VERSION = {v13},
  66. *}
  67. *
  68. *@article{franzina2010existence,
  69. * title={Existence and uniqueness for a p-laplacian nonlinear eigenvalue problem},
  70. * author={Franzina, Giovanni and Lamberti, Pier Domenico},
  71. * journal={Electron. J. Differential Equations},
  72. * volume={26},
  73. * number={10},
  74. * year={2010}
  75. *}
  76. *
  77. *
  78. *@article{belyaev,
  79. *author = {Belyaev, Alexander G. and Fayolle, Pierre-Alain},
  80. *title = {On Variational and PDE-Based Distance Function Approximations},
  81. *journal = {Computer Graphics Forum},
  82. *volume = {34},
  83. *number = {8},
  84. *issn = {1467-8659},
  85. *url = {http://dx.doi.org/10.1111/cgf.12611},
  86. *doi = {10.1111/cgf.12611},
  87. *pages = {104--118},
  88. *year = {2015},
  89. *}
  90. *
  91. *
  92. *
  93. * LANGAGE : GIBIANE-CAST3M
  94. * AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SEMT/LTA)
  95. * mél : stephane.gounand@cea.fr
  96. **********************************************************************
  97. * VERSION : v1, 02/05/2017, version initiale
  98. * HISTORIQUE : v1, 02/05/2017, création
  99. * HISTORIQUE :
  100. * HISTORIQUE :
  101. ************************************************************************
  102. 'OPTI' 'ECHO' 0 ;
  103. 'SAUT' 'PAGE' ;
  104. interact = faux ;
  105. graph = faux ;
  106. *
  107. * Procédure d'accélération de la convergence.
  108. *
  109. 'DEBP' ACCCVG ;
  110. 'ARGU' res*'CHPOINT' ;
  111. 'ARGU' tabacc*'TABLE' ;
  112. 'SI' ('NON' ('EXIS' tabacc 'iter')) ;
  113. tabacc . 'iter' = 0 ;
  114. tabacc . 'znacce' = 4 ;
  115. tabacc . 'itdep' = 3 ;
  116. tabacc . 'acfp1' = 'COPI' ('*' res 0.) ;
  117. tabacc . 'acfp2' = tabacc . 'acfp1' ;
  118. tabacc . 'acfp3' = tabacc . 'acfp1' ;
  119. tabacc . 'acfep1' = tabacc . 'acfp1' ;
  120. tabacc . 'acfep2' = tabacc . 'acfp1' ;
  121. tabacc . 'correc' = 0. ;
  122. tabacc . 'freap' = 0. ;
  123. 'FINS' ;
  124. it = '+' (tabacc . 'iter') 1 ;
  125. znacce = tabacc . 'znacce' ;
  126. itdep = tabacc . 'itdep' ;
  127. tabacc . 'iter' = it ;
  128. tabacc . 'correcp' = tabacc . 'correc' ;
  129. tabacc . 'correc' = 0. ;
  130. tabacc . 'acfp0' = 'ENLEVER' ('-' res (tabacc . 'freap')) 'FLX' ;
  131. tabacc . 'acfep0' = '-' (tabacc . 'acfp0') (tabacc . 'correcp') ;
  132. 'SI' ('ET' ('MULT' it znacce) ('>' it itdep)) ;
  133. 'MESS' 'acceleration' ;
  134. correc = 'ACT3' (tabacc . 'acfep2') (tabacc . 'acfep1')
  135. (tabacc . 'acfep0') (tabacc . 'acfp3')
  136. (tabacc . 'acfp2') (tabacc . 'acfp1')
  137. (tabacc . 'acfp0') ;
  138. res = '-' res correc ;
  139. 'FINS' ;
  140. 'SI' ('>' it 3) ;
  141. 'DETR' (tabacc . 'acfp3') ;
  142. 'DETR' (tabacc . 'acfep2') ;
  143. 'FINSI' ;
  144. tabacc . 'acfp3' = tabacc . 'acfp2 ' ;
  145. tabacc . 'acfp2' = tabacc . 'acfp1 ' ;
  146. tabacc . 'acfp1' = tabacc . 'acfp0 ' ;
  147. tabacc . 'acfep2' = tabacc . 'acfep1 ' ;
  148. tabacc . 'acfep1' = tabacc . 'acfep0 ' ;
  149. 'RESP' res tabacc ;
  150. 'FINP' ;
  151. *
  152. * Options globales
  153. *
  154. 'OPTI' 'DIME' 2 ;
  155. 'OPTI' 'ELEM' 'TRI3' ;
  156. *
  157. * Maillage
  158. *
  159. * nx : nombre de maille par côté
  160. * idom : type de domaine 0 (carré) ou 1 (quart de cercle)
  161. * Sur 1, il y a une solution analytique pour p>1
  162. nx = 40 ;
  163. idom = 1 ;
  164. * Physique
  165. * vp : exposant p du p-laplacien
  166. vp = 3. ;
  167. * Numérique
  168. * nitmax : nombre d'itérations max.
  169. * tol : tolérance pour la convergence de la boucle non-linéaire
  170. * omega : facteur de relaxation
  171. * iktan : type de matrice tangente :
  172. * 0 laplacien ; 1 point fixe ou ; 2 newton
  173. nitmax = 100 ;
  174. tol = 1.D-10 ;
  175. *omega = '/' 2. vp ;
  176. *iktan = 1 ;
  177. omega = 1. ;
  178. iktan = 2 ;
  179. lacc = 0 ; tabacc = 'TABL' ;
  180.  
  181. chpar = 'CHAI' 'p=' vp ' omega=' omega ' iktan=' iktan ' lacc=' lacc ;
  182.  
  183. 'MESS' '*************************' ;
  184. 'MESS' '* P laplacien ' chpar ;
  185. 'MESS' '*************************' ;
  186. * Points
  187. 'OPTI' 'DENS' ('/' 1. ('FLOT' nx)) ;
  188. pA = 0. 0. ;
  189. pB = 1. 0. ;
  190. pC = 1. 1. ;
  191. pD = 0. 1. ;
  192. *
  193. * Maillage carré
  194. *
  195. 'SI' ('EGA' idom 0) ;
  196. * Lignes
  197. lAB = 'DROI' nx pA pB ;
  198. lBC = 'DROI' nx pB pC ;
  199. lCD = 'DROI' nx pC pD ;
  200. lDA = 'DROI' nx pD pA ;
  201. * Surface
  202. *mt = 'DALL' lAB lBC lCD lDA ;
  203. mt = 'SURF' (lAB 'ET' lBC 'ET' lCD 'ET' lDA) ;
  204. * Blocages
  205. mailblo = 'CONT' mt ;
  206. 'FINS' ;
  207. *
  208. * Maillage quart de cercle
  209. *
  210. 'SI' ('EGA' idom 1) ;
  211. * Lignes
  212. lAB = 'DROI' pA pB ;
  213. lBD = 'CERC' pB pA pD ;
  214. lDA = 'DROI' pD pA ;
  215. * Surface
  216. *mt = 'DALL' lAB lBC lCD lDA ;
  217. mt = 'SURF' (lAB 'ET' lBD 'ET' lDA) ;
  218. mailblo = lBD ;
  219. 'FINS' ;
  220. *
  221. 'SI' graph ;
  222. 'TRAC' mt 'TITR' 'MAILLAGE' ;
  223. 'FINS' ;
  224. *
  225. * Modele et materiau initial
  226. *
  227. modt = 'MODE' mt 'THERMIQUE' 'ISOTROPE' ;
  228. modt2 = 'MODE' mt 'THERMIQUE' 'ANISOTROPE' ;
  229. matt = 'MATE' modt 'K' 1. ;
  230. *
  231. mlap0 = 'COND' modt matt ;
  232. * terme source
  233. tsou = 'SOUR' modt 1. mt ;
  234. * conditions aux limites
  235. mcli = 'BLOQ' 'T' mailblo ;
  236. * résolution
  237. mat = mlap0 'ET' mcli ;
  238. smb = tsou ;
  239. sol = 'RESO' mat smb ;
  240. solt = 'EXCO' 'T' sol 'T' ;
  241. *ch0 = 'Solution au pas 0' ;
  242. *'TRAC' solt mt 'TITR' ch0 ;
  243. ipas = 0 ; soli = sol ;
  244. lnres = 'PROG' ;
  245. cgtx = 'MOTS' 'T,X' ; cgty = 'MOTS' 'T,Y' ;
  246. cgt = cgtx 'ET' cgty ;
  247. ck11 = 'MOTS' 'K11' ; ck21 = 'MOTS' 'K21' ; ck22 = 'MOTS' 'K22' ;
  248. *
  249. * Boucle pour la non-linearite
  250. *
  251. 'REPE' bcl nitmax ;
  252. gsoli = 'GRAD' modt soli ;
  253. ngsoli = 'PSCA' gsoli gsoli cgt cgt ;
  254. * Nouveau matériau avec la nouvelle conductivité
  255. chk = ngsoli '**' ('/' ('-' vp 2.) 2.) ;
  256. matt = 'MATE' modt 'K' chk ;
  257. * Matrice et calcul du résidu
  258. mlap = 'COND' modt matt ;
  259. mat = mlap 'ET' mcli ;
  260. res = tsou '-' ('*' mat soli) ;
  261. *
  262. nres = 'MAXI' res 'ABS' 'AVEC' ('MOTS' 'Q') ;
  263. lnres = 'ET' lnres ('PROG' nres) ;
  264. 'SI' ('<' nres tol) ;
  265. 'MESS' ('CHAI' 'Pas i=' ipas ' max|Ri|=' nres ' < tol=' tol) ;
  266. 'QUIT' bcl ;
  267. 'SINO' ;
  268. 'MESS' ('CHAI' 'Pas i=' ipas ' max|Ri|=' nres) ;
  269. 'FINS' ;
  270. * Matrice tangente = rigidite au pas 0
  271. 'SI' ('EGA' iktan 0) ;
  272. mattan = mlap0 ;
  273. 'FINS' ;
  274. * Matrice tangente = Matrice de point fixe
  275. 'SI' ('EGA' iktan 1) ;
  276. mattan = mlap ;
  277. 'FINS' ;
  278. 'SI' ('EGA' iktan 2) ;
  279. * Matrice tangente = Matrice de Newton
  280. * 2 \eta' (|grad T|^2) (grad T \ptens grad T)
  281. * avec \eta'(z) = \frac{p-2}{2} z^{\frac{p-4}{2}}
  282. coeff = '*' (ngsoli '**' ('/' ('-' vp 4.) 2.))
  283. ('-' vp 2.) ;
  284. * On change le type de gsoli sinon il veut une syntaxe avec un modele
  285. gsoli = 'CHAN' 'TYPE' gsoli ' ' ;
  286. vk11 = '*' gsoli gsoli cgtx cgtx ck11 ;
  287. vk21 = '*' gsoli gsoli cgty cgtx ck21 ;
  288. vk22 = '*' gsoli gsoli cgty cgty ck22 ;
  289. matt2 = 'MATE' modt2 'DIRECTION' (1. 0.) 'PARALLELE'
  290. 'K11' ('*' vk11 coeff)
  291. 'K21' ('*' vk21 coeff)
  292. 'K22' ('*' vk22 coeff) ;
  293. matnew2 = 'COND' modt2 matt2 ;
  294. mattan = mlap 'ET' matnew2 ;
  295. 'FINS' ;
  296. mattot = ('*' mattan ('/' 1. omega)) 'ET' mcli ;
  297. smb = res ;
  298. 'SI' ('EGA' lacc 1) ; smb tabacc = ACCCVG smb tabacc ; 'FINS' ;
  299. dsol = 'RESO' mattot smb ;
  300. 'SI' ('EGA' lacc 1) ; tabacc . 'freap' = '*' mcli dsol ; 'FINS' ;
  301. solip = '+' soli dsol ;
  302. * solipt = 'EXCO' 'T' solip 'T' ;
  303. * chip = 'CHAI' 'Solution vp=' vp ' au pas ' ('+' ipas 1) ;
  304. * 'TRAC' solipt mt 'TITR' chip 'NCLK' ;
  305. soli = solip ;
  306. ipas = '+' ipas 1 ;
  307. 'FIN' bcl ;
  308. *
  309. * Courbe du résidu
  310. *
  311. npas = 'DIME' lnres ;
  312. 'SI' ('>' npas 1) ;
  313. lpas = 'PROG' 0. 'PAS' 1. 'NPAS' ('-' npas 1) ;
  314. llnres = '/' ('LOG' lnres) ('LOG' 10.) ;
  315. evl = 'EVOL' 'MANU' lpas llnres ;
  316. 'SI' graph ;
  317. titx = 'iter' ; tity = 'Log. |res|' ;
  318. tit = 'CHAI' tity '(' titx ')' ' ' chpar ;
  319. 'DESS' evl 'TITR' tit 'TITX' titx 'TITY' tity ;
  320. 'FINS' ;
  321. 'FINS' ;
  322. *
  323. * Tracé en élévation
  324. *
  325. 'SI' graph ;
  326. 'OPTI' 'DIME' 3 ;
  327. solit = 'EXCO' 'T' soli ;
  328. depz = 'NOMC' solit 'UZ' ;
  329. 'FORM' ('/' depz ('MAXI' depz 'ABS')) ;
  330. 'OPTI' 'OEIL' (20. 30. 10.) ;
  331. 'OPTI' 'ISOV' 'SULI' ;
  332. chi = 'CHAI' 'Solution vp=' vp ' au pas ' ipas ;
  333. 'TRAC' solit mt ('CONT' mt) 30 'TITR' chi ;
  334. 'OPTI' 'DIME' 2 ;
  335. 'FINS' ;
  336. *
  337. * Comparaison à la solution analytique
  338. *
  339. 'SI' ('EGA' idom 1) ;
  340. x y = 'COOR' mt ;
  341. r2 = (x '*' x) '+' (y '*' y) ;
  342. expo = vp '/' ('-' vp 1.) ;
  343. vd = 'VALE' 'DIME' ;
  344. Const = '/' (vd '**' ('/' -1 ('-' vp 1.)))
  345. expo ;
  346. solana = '*' ('-' 1. ('**' r2 ('/' expo 2.)))
  347. const ;
  348. solnum = 'EXCO' 'T' soli ;
  349. dsol = '-' solana solnum ;
  350. 'SI' graph ;
  351. 'TRAC' dsol mt ('CONT' mt) 30
  352. 'TITR' 'Difference analytique-numerique' ;
  353. 'FINS' ;
  354. mdsol = 'MAXI' dsol 'ABS' ;
  355. ch = 'CHAI' 'Maxi. |analytique - numerique|=' mdsol ' ' chpar;
  356. 'MESS' ch ;
  357. 'FINS' ;
  358. *
  359. * Tests :
  360. * 1) La méthode de Newton converge en 7-8 itérations pour p=3 (cf. Saramito
  361. * p.137)
  362. * 2) L'erreur par rapport à la solution analytique vaut environ 4.e-4
  363. * avec 40 mailles
  364. *
  365. lok = vrai ;
  366. lok = 'ET' lok (npas '&lt;EG' 8) ;
  367. 'SI' ('NON' lok) ;
  368. ch = 'CHAI' 'Newton : unexpectedly slow convergence' ;
  369. 'ERRE' ch ;
  370. 'FINS' ;
  371. lok = 'ET' lok (mdsol '<' 5.e-4) ;
  372. 'SI' ('NON' lok) ;
  373. ch = 'CHAI' 'Comparison with analytical solution : bad results' ;
  374. 'ERRE' ch ;
  375. 'FINS' ;
  376. 'SI' ('NON' lok) ;
  377. 'ERREUR' 5 ;
  378. 'SINON' ;
  379. 'MESSAGE' ('CHAINE' 'Test ok !') ;
  380. 'FINSI' ;
  381. *
  382. 'SI' interact ;
  383. 'OPTION' 'DONN' 5 'ECHO' 1 ;
  384. 'FINSI' ;
  385. *
  386. * End of dgibi file THER9
  387. *
  388. 'FIN' ;
  389.  
  390.  
  391.  
  392.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales