1 : C YLAP1C SOURCE MAGN 09/07/30 21:22:17 6458
2 : SUBROUTINE YLAP1C(ICOGRV,MPROC,MPVITC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,LCLIMV,KRVIMP,LCLITO,KRTOIM,
5 : $ NOMINC,
6 : $ MU,
7 : $ IJACO,
8 : $ IMPR,IRET)
9 : IMPLICIT INTEGER(I-N)
10 : IMPLICIT REAL*8 (A-H,O-Z)
11 : C***********************************************************************
12 : C NOM : YLAP1C
13 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
14 : C VF 2D.
15 : C Ici, on ne calcule que les contributions 'simples'
16 : C à la matrice jacobienne faisant intervenir les
17 : C coefficients pour le calcul des gradients de vitesse
18 : C (ICOGRV).
19 : C (contributions à (d Res_{\rho u} / d var) et (d
20 : C Res_{\rho v} / d var)
21 : C var prenant successivement les valeurs :
22 : C \rho, \rho u, \rho v, \rho e_t )
23 : C
24 : C
25 : C NOTE PERSO : On pourrait programmer ça plus lisiblement en stockant
26 : C les contributions dans un tableau de pointeurs (2
27 : C indices, c'est possible ?) et en effectuant des produits
28 : C matriciels (coeff. x matrices de dérivées).
29 : C On n'y coupera pas si on veut regrouper 2D et 3D...
30 : C On ne va pas le faire.
31 : C
32 : C
33 : C LANGAGE : ESOPE
34 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
35 : C mél : gounand@semt2.smts.cea.fr
36 : C***********************************************************************
37 : C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non
38 : C standard) à un objet de type MATRIK.
39 : C APPELE PAR : YLAP1A : Calcul de la matrice jacobienne du
40 : C résidu du laplacien VF 2D.
41 : C***********************************************************************
42 : C ENTREES : ICOGRV (type MCHELM) : coefficients pour le
43 : C calcul du gradient de la vitesse aux
44 : C interfaces.
45 : C MPROC (type MPOVAL) : masse volumique par
46 : C élément.
47 : C MPVITC (type MPOVAL) : vitesse par élément.
48 : C MPVOLU (type MPOVAL) : volume des éléments.
49 : C MPNORM (type MPOVAL) : normale aux faces.
50 : C MPSURF (type MPOVAL) : surface des faces.
51 : C MELEFL (type MELEME) : connectivités face-(centre
52 : C gauche, centre droit).
53 : C KRFACE (type MLENTI) : tableau de repérage dans
54 : C le maillage des faces des éléments.
55 : C KRCENT (type MLENTI) : tableau de repérage dans
56 : C le maillage des centres des éléments.
57 : C LCLIMV (type logique) : .TRUE. => CL de Dirichlet
58 : C sur la vitesse.
59 : C KRVIMP (type MLENTI) : tableau de repérage dans
60 : C maillage des CL de Dirichlet sur la
61 : C vitesse.
62 : C LCLITO (type logique) : .TRUE. => CL de Dirichlet
63 : C sur le tenseur des contraintes.
64 : C KRTOIM (type MLENTI) : tableau de repérage dans
65 : C maillage des CL de Dirichlet sur le tenseur des
66 : C contraintes
67 : C NOMINC (type MLMOTS) : noms des inconnues.
68 : C MU (type réel) : viscosité dynamique (SI).
69 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
70 : C résidu du laplacien VF 2D.
71 : C SORTIES : -
72 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
73 : C***********************************************************************
74 : C VERSION : v1, 08/08/2001, version initiale
75 : C HISTORIQUE : v1, 08/08/2001, création
76 : C HISTORIQUE :
77 : C HISTORIQUE :
78 : C***********************************************************************
79 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
80 : C en cas de modification de ce sous-programme afin de faciliter
81 : C la maintenance !
82 : C***********************************************************************
83 : -INC CCOPTIO
84 : -INC SMCHPOI
85 : POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL
86 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
87 : -INC SMCHAML
88 : POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML
89 : POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL
90 : -INC SMELEME
91 : POINTEUR MELEFL.MELEME
92 : POINTEUR MCOGRV.MELEME
93 : -INC SMLENTI
94 : POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI
95 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
96 : -INC SMLMOTS
97 : POINTEUR NOMINC.MLMOTS
98 : -INC SMMATRIK
99 : POINTEUR IJACO.MATRIK
100 : *
101 : * Objet matrice élémentaire simplifié
102 : *
103 : SEGMENT GMATSI
104 : INTEGER POIPR1(NPP1,NEL1)
105 : INTEGER POIDU1(1,NEL1)
106 : INTEGER POIPR2(NPP2,NEL2)
107 : INTEGER POIDU2(2,NEL2)
108 : POINTEUR LMATSI(0).MATSIM
109 : ENDSEGMENT
110 : * Contributions simples de la part du gradient de
111 : * vitesse (CTSGRV)
112 : POINTEUR CTSGRV.GMATSI
113 : SEGMENT MATSIM
114 : CHARACTER*8 NOMPRI,NOMDUA
115 : REAL*8 VALMA1(1,NPP1,NEL1)
116 : REAL*8 VALMA2(2,NPP2,NEL2)
117 : ENDSEGMENT
118 : POINTEUR ROURHO.MATSIM
119 : POINTEUR ROVRHO.MATSIM
120 : POINTEUR ROUROU.MATSIM
121 : POINTEUR ROUROV.MATSIM
122 : POINTEUR ROVROU.MATSIM
123 : POINTEUR ROVROV.MATSIM
124 : *
125 : REAL*8 MU
126 : *
127 : INTEGER IMPR,IRET
128 : *
129 : LOGICAL LCLIMV,LCLITO
130 : LOGICAL LMUR
131 : LOGICAL LCTRB1,LCTRB2
132 : *
133 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
134 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
135 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
136 : INTEGER NLCENP,NLCEND,NLFACE,NLCLTO,NLCLV
137 : INTEGER NPTEL
138 : *
139 : REAL*8 ALPHAX,ALPHAY,CNX,CNY
140 : REAL*8 SIGNOR,SURFFA,VOLUEL
141 : REAL*8 RHOP,UP,VP
142 : REAL*8 FACTOR
143 : REAL*8 DUDRHO,DUDROU,DVDRHO,DVDROV
144 : REAL*8 BROUDU,BROUDV,BROVDU,BROVDV
145 : C
146 : *
147 : * Executable statements
148 : *
149 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap1c.eso'
150 : * On calcule les contributions à (d Res_{\rho u} / d var) et (d
151 : * Res_{\rho v} / d var) ; var prenant successivement les valeurs :
152 : * \rho, \rho u, \rho v, \rho e_t
153 : * Note :
154 : * - (d Res_{\rho u} / d \rho e_t) = 0
155 : * - (d Res_{\rho v} / d \rho e_t) = 0
156 : * On dérive les termes : \tens{\tau} \prod \vect{n}
157 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
158 : * ROURHO, ROUROU, ROUROV, ROVRHO, ROVROU, ROVROV
159 : IF (LCLIMV) THEN
160 : SEGACT KRVIMP
161 : ENDIF
162 : IF (LCLITO) THEN
163 : SEGACT KRTOIM
164 : ENDIF
165 : SEGACT NOMINC
166 : SEGACT KRCENT
167 : SEGACT KRFACE
168 : SEGACT MELEFL
169 : SEGACT MPSURF
170 : SEGACT MPNORM
171 : SEGACT MPVOLU
172 : SEGACT MPROC
173 : SEGACT MPVITC
174 : SEGACT ICOGRV
175 : NSOUCH=ICOGRV.IMACHE(/1)
176 : DO 1 ISOUCH=1,NSOUCH
177 : MCOGRV=ICOGRV.IMACHE(ISOUCH)
178 : JCOGRV=ICOGRV.ICHAML(ISOUCH)
179 : SEGACT JCOGRV
180 : KDUNDX=JCOGRV.IELVAL(1)
181 : KDUNDY=JCOGRV.IELVAL(2)
182 : SEGDES JCOGRV
183 : SEGACT KDUNDX
184 : SEGACT KDUNDY
185 : SEGACT MCOGRV
186 : NELEM=MCOGRV.NUM(/2)
187 : NPTEL=MCOGRV.NUM(/1)
188 : NPP1=NPTEL-1
189 : NPP2=NPTEL-1
190 : NEL1=NELEM
191 : NEL2=NELEM
192 : IEL1=1
193 : IEL2=1
194 : SEGINI ROURHO
195 : SEGINI ROVRHO
196 : SEGINI ROUROU
197 : SEGINI ROVROU
198 : SEGINI ROUROV
199 : SEGINI ROVROV
200 : SEGINI CTSGRV
201 : ROURHO.NOMPRI(1:4)=NOMINC.MOTS(1)
202 : ROURHO.NOMPRI(5:8)=' '
203 : ROURHO.NOMDUA(1:4)=NOMINC.MOTS(2)
204 : ROURHO.NOMDUA(5:8)=' '
205 : ROVRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
206 : ROVRHO.NOMPRI(5:8)=' '
207 : ROVRHO.NOMDUA(1:4)=NOMINC.MOTS(3)
208 : ROVRHO.NOMDUA(5:8)=' '
209 : ROUROU.NOMPRI(1:4)=NOMINC.MOTS(2)
210 : ROUROU.NOMPRI(5:8)=' '
211 : ROUROU.NOMDUA(1:4)=NOMINC.MOTS(2)
212 : ROUROU.NOMDUA(5:8)=' '
213 : ROVROU.NOMPRI(1:4)=NOMINC.MOTS(2)
214 : ROVROU.NOMPRI(5:8)=' '
215 : ROVROU.NOMDUA(1:4)=NOMINC.MOTS(3)
216 : ROVROU.NOMDUA(5:8)=' '
217 : ROUROV.NOMPRI(1:4)=NOMINC.MOTS(3)
218 : ROUROV.NOMPRI(5:8)=' '
219 : ROUROV.NOMDUA(1:4)=NOMINC.MOTS(2)
220 : ROUROV.NOMDUA(5:8)=' '
221 : ROVROV.NOMPRI(1:4)=NOMINC.MOTS(3)
222 : ROVROV.NOMPRI(5:8)=' '
223 : ROVROV.NOMDUA(1:4)=NOMINC.MOTS(3)
224 : ROVROV.NOMDUA(5:8)=' '
225 : DO 12 IELEM=1,NELEM
226 : * Le premier point du support de ICOGRV est un point FACE
227 : NGFACE=MCOGRV.NUM(1,IELEM)
228 : NLFACE=KRFACE.LECT(NGFACE)
229 : IF (NLFACE.EQ.0) THEN
230 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
231 : GOTO 9999
232 : ENDIF
233 : * On calcule la contribution à la matrice jacobienne IJACO de la face
234 : * NGFAC (points duaux : centres à gauche et à droite de la face)
235 : * (points primaux : une partie (bicoz conditions aux limites)
236 : * de ceux du stencil pour le calcul du gradient
237 : * à la face, ils doivent être des points centres)
238 : * Si le tenseur des contraintes sur la face est imposé par les
239 : * conditions aux limites, la contribution de la face à IJACO est
240 : * nulle.
241 : LCTRB1=.TRUE.
242 : IF (LCLITO) THEN
243 : NLCLTO=KRTOIM.LECT(NGFACE)
244 : IF (NLCLTO.NE.0) THEN
245 : LCTRB1=.FALSE.
246 : ENDIF
247 : ENDIF
248 : IF (LCTRB1) THEN
249 : NGCGAU=MELEFL.NUM(1,NLFACE)
250 : NGCDRO=MELEFL.NUM(3,NLFACE)
251 : LMUR=(NGCGAU.EQ.NGCDRO)
252 : * On distingue le cas où la face est un bord du maillage (mur)
253 : * du cas où la face est interne au maillage
254 : IF (.NOT.LMUR) THEN
255 : NPD=2
256 : ELSE
257 : NPD=1
258 : ENDIF
259 : NPP=NPTEL-1
260 : * IPD=1 : point à gauche du point NGFACE
261 : * IPD=2 : point à droite du point NGFACE
262 : DO 122 IPD=1,NPD
263 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
264 : IF (.NOT.LMUR) THEN
265 : CTSGRV.POIDU2(IPD,IEL2)=NPDUAL
266 : ELSE
267 : CTSGRV.POIDU1(IPD,IEL1)=NPDUAL
268 : ENDIF
269 : NLCEND=KRCENT.LECT(NPDUAL)
270 : IF (NLCEND.EQ.0) THEN
271 : WRITE(IOIMP,*) 'Erreur grave n°1'
272 : GOTO 9999
273 : ENDIF
274 : DO 124 IPP=1,NPP
275 : NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
276 : C
277 : C* Modif AB: elimination du test condition limite
278 : C
279 : NLCENP=KRCENT.LECT(NPPRIM)
280 : IF(NLCENP .EQ. 0) THEN
281 : * Lorsque une contribution est nulle, on fixe artificiellement le
282 : * point primal égal au point dual.
283 : IF (.NOT.LMUR) THEN
284 : CTSGRV.POIPR2(IPP,IEL2)=NPDUAL
285 : ROURHO.VALMA2(IPD,IPP,IEL2)=0.D0
286 : ROVRHO.VALMA2(IPD,IPP,IEL2)=0.D0
287 : ROUROU.VALMA2(IPD,IPP,IEL2)=0.D0
288 : ROVROU.VALMA2(IPD,IPP,IEL2)=0.D0
289 : ROUROV.VALMA2(IPD,IPP,IEL2)=0.D0
290 : ROVROV.VALMA2(IPD,IPP,IEL2)=0.D0
291 : ELSE
292 : CTSGRV.POIPR1(IPP,IEL1)=NPDUAL
293 : ROURHO.VALMA1(IPD,IPP,IEL1)=0.D0
294 : ROVRHO.VALMA1(IPD,IPP,IEL1)=0.D0
295 : ROUROU.VALMA1(IPD,IPP,IEL1)=0.D0
296 : ROVROU.VALMA1(IPD,IPP,IEL1)=0.D0
297 : ROUROV.VALMA1(IPD,IPP,IEL1)=0.D0
298 : ROVROV.VALMA1(IPD,IPP,IEL1)=0.D0
299 : ENDIF
300 : ELSE
301 : * Les contributions valent :
302 : * 1. (d Res_{\rho u})_d / (d var)_p =
303 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \mu
304 : * * [ [ ( 4/3 (n_x * \alpha_x) + (n_y * \alpha_y)) *
305 : * ((du)_p / (d var)_p)]
306 : * + [ (-2/3 (n_x * \alpha_y) + (n_y * \alpha_x)) *
307 : * ((dv)_p / (d var)_p)]
308 : * ]
309 : * 2. (d Res_{\rho v})_d / (d var)_p =
310 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \mu
311 : * * [ [ (-2/3 (n_y * \alpha_x) + (n_x * \alpha_y)) *
312 : * ((du)_p / (d var)_p)]
313 : * + [ ( 4/3 (n_y * \alpha_y) + (n_x * \alpha_x)) *
314 : * ((dv)_p / (d var)_p)]
315 : * ]
316 : *
317 : * avec :
318 : * (du)_p / (d \rho)_p = - u / \rho_p
319 : * (du)_p / (d \rho u)_p = 1 / \rho_p
320 : * (du)_p / (d \rho v)_p = 0
321 : * (dv)_p / (d \rho)_p = - v / \rho_p
322 : * (dv)_p / (d \rho u)_p = 0
323 : * (dv)_p / (d \rho v)_p = 1 / \rho_p
324 : *
325 : IF (NLCENP.EQ.0) THEN
326 : WRITE(IOIMP,*) 'Erreur grave n°2'
327 : GOTO 9999
328 : ENDIF
329 : * normale sortante pour IPD=1, rentrante pour IPD=2
330 : SIGNOR=(-1.D0)**(IPD+1)
331 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
332 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
333 : CNX =MPNORM.VPOCHA(NLFACE,1)
334 : CNY =MPNORM.VPOCHA(NLFACE,2)
335 : ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
336 : ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
337 : RHOP =MPROC.VPOCHA(NLCENP,1)
338 : UP =MPVITC.VPOCHA(NLCENP,1)
339 : VP =MPVITC.VPOCHA(NLCENP,2)
340 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*MU
341 : BROUDU=(( 4.D0/3.D0)*(CNX*ALPHAX))
342 : $ + (CNY*ALPHAY)
343 : BROUDV=((-2.D0/3.D0)*(CNX*ALPHAY))
344 : $ + (CNY*ALPHAX)
345 : BROVDU=((-2.D0/3.D0)*(CNY*ALPHAX))
346 : $ + (CNX*ALPHAY)
347 : BROVDV=(( 4.D0/3.D0)*(CNY*ALPHAY))
348 : $ + (CNX*ALPHAX)
349 : DUDRHO=-UP /RHOP
350 : DUDROU=1.D0/RHOP
351 : DVDRHO=-VP /RHOP
352 : DVDROV=1.D0/RHOP
353 : IF (.NOT.LMUR) THEN
354 : CTSGRV.POIPR2(IPP,IEL2)=NPPRIM
355 : ROURHO.VALMA2(IPD,IPP,IEL2)=
356 : $ FACTOR*((BROUDU*DUDRHO)+(BROUDV*DVDRHO))
357 : ROVRHO.VALMA2(IPD,IPP,IEL2)=
358 : $ FACTOR*((BROVDU*DUDRHO)+(BROVDV*DVDRHO))
359 : ROUROU.VALMA2(IPD,IPP,IEL2)=
360 : $ FACTOR* (BROUDU*DUDROU)
361 : ROVROU.VALMA2(IPD,IPP,IEL2)=
362 : $ FACTOR* (BROVDU*DUDROU)
363 : ROUROV.VALMA2(IPD,IPP,IEL2)=
364 : $ FACTOR* (BROUDV*DVDROV)
365 : ROVROV.VALMA2(IPD,IPP,IEL2)=
366 : $ FACTOR* (BROVDV*DVDROV)
367 : ELSE
368 : CTSGRV.POIPR1(IPP,IEL1)=NPPRIM
369 : ROURHO.VALMA1(IPD,IPP,IEL1)=
370 : $ FACTOR*((BROUDU*DUDRHO)+(BROUDV*DVDRHO))
371 : ROVRHO.VALMA1(IPD,IPP,IEL1)=
372 : $ FACTOR*((BROVDU*DUDRHO)+(BROVDV*DVDRHO))
373 : ROUROU.VALMA1(IPD,IPP,IEL1)=
374 : $ FACTOR* (BROUDU*DUDROU)
375 : ROVROU.VALMA1(IPD,IPP,IEL1)=
376 : $ FACTOR* (BROVDU*DUDROU)
377 : ROUROV.VALMA1(IPD,IPP,IEL1)=
378 : $ FACTOR* (BROUDV*DVDROV)
379 : ROVROV.VALMA1(IPD,IPP,IEL1)=
380 : $ FACTOR* (BROVDV*DVDROV)
381 : ENDIF
382 : ENDIF
383 : 124 CONTINUE
384 : 122 CONTINUE
385 : IF (.NOT.LMUR) THEN
386 : IEL2=IEL2+1
387 : ELSE
388 : IEL1=IEL1+1
389 : ENDIF
390 : ENDIF
391 : 12 CONTINUE
392 : NPP1=NPTEL-1
393 : NPP2=NPTEL-1
394 : NEL1=IEL1-1
395 : NEL2=IEL2-1
396 : SEGADJ ROURHO
397 : SEGADJ ROVRHO
398 : SEGADJ ROUROU
399 : SEGADJ ROVROU
400 : SEGADJ ROUROV
401 : SEGADJ ROVROV
402 : SEGADJ CTSGRV
403 : CTSGRV.LMATSI(**)=ROURHO
404 : CTSGRV.LMATSI(**)=ROVRHO
405 : CTSGRV.LMATSI(**)=ROUROU
406 : CTSGRV.LMATSI(**)=ROVROU
407 : CTSGRV.LMATSI(**)=ROUROV
408 : CTSGRV.LMATSI(**)=ROVROV
409 : * On accumule les matrices résultantes dans IJACO
409 : CALL AJMTK(CTSGRV,IJACO,IMPR,IRET)
410 : IF (IRET.NE.0) GOTO 9999
411 : SEGSUP ROURHO
412 : SEGSUP ROVRHO
413 : SEGSUP ROUROU
414 : SEGSUP ROVROU
415 : SEGSUP ROUROV
416 : SEGSUP ROVROV
417 : SEGSUP CTSGRV
418 : *
419 : SEGDES MCOGRV
420 : SEGDES KDUNDY
421 : SEGDES KDUNDX
422 : 1 CONTINUE
423 : SEGDES ICOGRV
424 : SEGDES MPVITC
425 : SEGDES MPROC
426 : SEGDES MPVOLU
427 : SEGDES MPNORM
428 : SEGDES MPSURF
429 : SEGDES MELEFL
430 : SEGDES KRFACE
431 : SEGDES KRCENT
432 : SEGDES NOMINC
433 : IF (LCLITO) THEN
434 : SEGDES KRTOIM
435 : ENDIF
436 : IF (LCLIMV) THEN
437 : SEGDES KRVIMP
438 : ENDIF
439 : *
440 : * Normal termination
441 : *
442 : IRET=0
443 : RETURN
444 : *
445 : * Format handling
446 : *
447 : *
448 : * Error handling
449 : *
450 : 9999 CONTINUE
451 : IRET=1
452 : WRITE(IOIMP,*) 'An error was detected in subroutine ylap1c'
453 : RETURN
454 : *
455 : * End of subroutine YLAP1C
456 : *
457 : END
458 :
459 :
460 :
461 :
462 :
463 :
464 :
465 :
© Cast3M 2003 - All rights reserved.
Disclaimer