1 : C ZLAP2B SOURCE MAGN 09/07/30 21:23:43 6458
2 : SUBROUTINE ZLAP2B(PROPHY,PROPH2,MPROC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,LCLIMR,KRRIMP,MPRIMP,
5 : $ NOMINC,
6 : $ IJACO,
7 : $ IMPR,IRET)
8 : IMPLICIT INTEGER(I-N)
9 : IMPLICIT REAL*8 (A-H,O-Z)
10 : C***********************************************************************
11 : C NOM : ZLAP2B
12 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
13 : C VF 3D (termes multi-espèces).
14 : C Ici, on ne calcule que les contributions à la matrice
15 : C jacobienne faisant intervenir les coefficients pour le
16 : C calcul des gradients de Yk
17 : C (contributions à d Res_{\rho Yk} / d var
18 : C var prenant successivement les valeurs :
19 : C \rho, \rho Yk)
20 : C
21 : C
22 : C
23 : C
24 : C LANGAGE : ESOPE
25 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
26 : C mél : gounand@semt2.smts.cea.fr
27 : C***********************************************************************
28 : C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non
29 : C standard) à un objet de type MATRIK.
30 : C APPELE PAR : ZLAP2A : Calcul de la matrice jacobienne du
31 : C résidu du laplacien VF 3D.
32 : C***********************************************************************
33 : C ENTREES : PROPHY (type PROPHY) : propriétés des espèces
34 : C PROPH2 (type PROPH2) : précond. de PROPHY
35 : C MPROC (type MPOVAL) : masse volumique par
36 : C élément.
37 : C MPVOLU (type MPOVAL) : volume des éléments.
38 : C MPNORM (type MPOVAL) : normale aux faces.
39 : C MPSURF (type MPOVAL) : surface des faces.
40 : C MELEFL (type MELEME) : connectivités face-(centre
41 : C gauche, centre droit).
42 : C KRFACE (type MLENTI) : tableau de repérage dans
43 : C le maillage des faces des éléments.
44 : C KRCENT (type MLENTI) : tableau de repérage dans
45 : C le maillage des centres des éléments.
46 : C LCLIMR (type logique) : .TRUE. => CL de Dirichlet
47 : C sur la densité.
48 : C KRRIMP (type MLENTI) : tableau de repérage dans
49 : C maillage des CL de Dirichlet sur la densité.
50 : C MPRIMP (type MPOVAL) : valeurs des CL de
51 : C Dirichlet sur la densité.
52 : C NOMINC (type MLMOTS) : noms des inconnues.
53 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
54 : C résidu du laplacien VF 3D.
55 : C SORTIES : -
56 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
57 : C***********************************************************************
58 : C VERSION : v1, 08/03/2002, version initiale
59 : C HISTORIQUE : v1, 08/03/2002, création
60 : C HISTORIQUE :
61 : C HISTORIQUE :
62 : C***********************************************************************
63 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
64 : C en cas de modification de ce sous-programme afin de faciliter
65 : C la maintenance !
66 : C***********************************************************************
67 : -INC CCOPTIO
68 : -INC SMCOORD
69 : -INC SMCHPOI
70 : POINTEUR MPDKC.MPOVAL
71 : POINTEUR MPROC.MPOVAL,MPYKC.MPOVAL
72 : POINTEUR MPRIMP.MPOVAL
73 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
74 : -INC SMCHAML
75 : POINTEUR ICOGRY.MCHELM,JCOGRY.MCHAML
76 : POINTEUR KDYKDX.MELVAL,KDYKDY.MELVAL,KDYKDZ.MELVAL
77 : -INC SMELEME
78 : POINTEUR MELEFL.MELEME
79 : POINTEUR MCOGRY.MELEME
80 : -INC SMLENTI
81 : POINTEUR KRYIMP.MLENTI,KRRIMP.MLENTI
82 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
83 : -INC SMLMOTS
84 : POINTEUR NOMINC.MLMOTS
85 : -INC SMMATRIK
86 : POINTEUR IJACO.MATRIK
87 : INTEGER NESP,IESP
88 : SEGMENT PROPHY
89 : CHARACTER*4 NOMESP(NESP+1)
90 : REAL*8 CV(NESP+1)
91 : REAL*8 R(NESP+1)
92 : REAL*8 H0K(NESP+1)
93 : POINTEUR CDIFF(NESP+1).MCHPOI
94 : POINTEUR YK(NESP+1).MCHPOI
95 : POINTEUR GRADYK(NESP+1).MCHPOI
96 : POINTEUR CGRYK(NESP+1).MCHELM
97 : POINTEUR CLYK(NESP+1).MCHPOI
98 : ENDSEGMENT
99 : SEGMENT PROPH2
100 : POINTEUR MPDIFF(NESP+1).MPOVAL
101 : POINTEUR MPVALY(NESP+1).MPOVAL
102 : POINTEUR MPGRAD(NESP+1).MPOVAL
103 : LOGICAL LCLIM(NESP+1)
104 : POINTEUR KRCLIM(NESP+1).MLENTI
105 : ENDSEGMENT
106 : *
107 : * Objet matrice élémentaire simplifié
108 : *
109 : SEGMENT GMATSI
110 : INTEGER POIPR1(NPP1,NEL1)
111 : INTEGER POIDU1(1,NEL1)
112 : INTEGER POIPR2(NPP2,NEL2)
113 : INTEGER POIDU2(2,NEL2)
114 : POINTEUR LMATSI(0).MATSIM
115 : ENDSEGMENT
116 : * Contributions de la part du gradient de YK (CTGRY)
117 : POINTEUR CTGRY.GMATSI
118 : SEGMENT MATSIM
119 : CHARACTER*8 NOMPRI,NOMDUA
120 : REAL*8 VALMA1(1,NPP1,NEL1)
121 : REAL*8 VALMA2(2,NPP2,NEL2)
122 : ENDSEGMENT
123 : POINTEUR RYKRHO.MATSIM
124 : POINTEUR RYKRYK.MATSIM
125 : *
126 : REAL*8 DKG,DKD,DKF,RHOG,RHOD,RHOF
127 : *
128 : INTEGER IMPR,IRET
129 : *
130 : LOGICAL LCLIMY,LCLIMR
131 : LOGICAL LMUR
132 : LOGICAL LCTRB2
133 : *
134 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
135 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
136 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
137 : INTEGER NLCENP,NLCEND,NLFACE,NLCLY,NLFRI
138 : INTEGER NPTEL
139 : *
140 : REAL*8 BETAX,BETAY,BETAZ,CNX,CNY,CNZ
141 : REAL*8 SIGNOR,SURFFA,VOLUEL
142 : REAL*8 RHOP,YKP
143 : REAL*8 FACTOR
144 : REAL*8 DYKDRO,DYKDRY
145 : *
146 : INTEGER ICOORX,NLCGAU,NLCDRO
147 : REAL*8 XF,YF,ZF,XG,YG,ZG,XFMXG,YFMYG,ZFMZG,DRG
148 : & ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH
149 : C
150 : *
151 : * Executable statements
152 : *
153 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans zlap2b.eso'
154 : SEGACT PROPHY
155 : NESP=PROPHY.CV(/1)-1
156 : SEGACT PROPH2
157 : DO IESP=1,NESP
158 : IF (LCLIMR) THEN
159 : SEGACT KRRIMP
160 : SEGACT MPRIMP
161 : ENDIF
162 : LCLIMY=PROPH2.LCLIM(IESP)
163 : IF (LCLIMY) THEN
164 : KRYIMP=PROPH2.KRCLIM(IESP)
165 : SEGACT KRYIMP
166 : ENDIF
167 : SEGACT NOMINC
168 : SEGACT KRCENT
169 : SEGACT KRFACE
170 : SEGACT MELEFL
171 : SEGACT MPSURF
172 : SEGACT MPNORM
173 : SEGACT MPVOLU
174 : MPDKC=PROPH2.MPDIFF(IESP)
175 : SEGACT MPDKC
176 : SEGACT MPROC
177 : MPYKC=PROPH2.MPVALY(IESP)
178 : SEGACT MPYKC
179 : ICOGRY=PROPHY.CGRYK(IESP)
180 : SEGACT ICOGRY
181 : NSOUCH=ICOGRY.IMACHE(/1)
182 : DO 1 ISOUCH=1,NSOUCH
183 : MCOGRY=ICOGRY.IMACHE(ISOUCH)
184 : JCOGRY=ICOGRY.ICHAML(ISOUCH)
185 : SEGACT JCOGRY
186 : KDYKDX=JCOGRY.IELVAL(1)
187 : KDYKDY=JCOGRY.IELVAL(2)
188 : KDYKDZ=JCOGRY.IELVAL(3)
189 : SEGDES JCOGRY
190 : SEGACT KDYKDX
191 : SEGACT KDYKDY
192 : SEGACT KDYKDZ
193 : SEGACT MCOGRY
194 : NELEM=MCOGRY.NUM(/2)
195 : NPTEL=MCOGRY.NUM(/1)
196 : NPP1=NPTEL-1
197 : NPP2=NPTEL-1
198 : NEL1=NELEM
199 : NEL2=NELEM
200 : IEL1=1
201 : IEL2=1
202 : SEGINI RYKRHO
203 : SEGINI RYKRYK
204 : SEGINI CTGRY
205 : RYKRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
206 : RYKRHO.NOMPRI(5:8)=' '
207 : RYKRHO.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
208 : RYKRHO.NOMDUA(5:8)=' '
209 : RYKRYK.NOMPRI(1:4)=NOMINC.MOTS(IDIM+2+IESP)
210 : RYKRYK.NOMPRI(5:8)=' '
211 : RYKRYK.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
212 : RYKRYK.NOMDUA(5:8)=' '
213 : DO 12 IELEM=1,NELEM
214 : * Le premier point du support de ICOGRY est un point FACE
215 : NGFACE=MCOGRY.NUM(1,IELEM)
216 : NLFACE=KRFACE.LECT(NGFACE)
217 : IF (NLFACE.EQ.0) THEN
218 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
219 : GOTO 9999
220 : ENDIF
221 : * On calcule la contribution à la matrice jacobienne IJACO de la face
222 : * NGFAC (points duaux : centres à gauche et à droite de la face)
223 : * (points primaux : une partie (bicoz conditions aux limites)
224 : * de ceux du stencil pour le calcul du gradient
225 : * à la face, ils doivent être des points centres)
226 : * Si le flux de chaleur sur la face est imposé par les conditions
227 : * aux limites, la contribution de la face à IJACO est nulle.
228 : NGCGAU=MELEFL.NUM(1,NLFACE)
229 : NGCDRO=MELEFL.NUM(3,NLFACE)
230 : NLCGAU=KRCENT.LECT(NGCGAU)
231 : NLCDRO=KRCENT.LECT(NGCDRO)
232 : LMUR=(NGCGAU.EQ.NGCDRO)
233 : * On distingue le cas où la face est un bord du maillage (mur)
234 : * du cas où la face est interne au maillage
235 : IF (.NOT.LMUR) THEN
236 : NPD=2
237 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
238 : XF = MCOORD.XCOOR(ICOORX)
239 : YF = MCOORD.XCOOR(ICOORX+1)
240 : ZF = MCOORD.XCOOR(ICOORX+2)
241 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
242 : XG = MCOORD.XCOOR(ICOORX)
243 : YG = MCOORD.XCOOR(ICOORX+1)
244 : ZG = MCOORD.XCOOR(ICOORX+2)
245 : XFMXG = XF - XG
246 : YFMYG = YF - YG
247 : ZFMZG = ZF - ZG
248 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
249 : ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
250 : XD = MCOORD.XCOOR(ICOORX)
251 : YD = MCOORD.XCOOR(ICOORX+1)
252 : ZD = MCOORD.XCOOR(ICOORX+2)
253 : XFMXD = XF - XD
254 : YFMYD = YF - YD
255 : ZFMZD = ZF - ZD
256 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
257 : ALPHA=DRG/(DRG+DRD)
258 : UMALPH= 1.0D0 - ALPHA
259 : ELSE
260 : NPD=1
261 : ALPHA=0.0D0
262 : UMALPH=1.0D0
263 : ENDIF
264 : IF (LCLIMR) THEN
265 : NLFRI=KRRIMP.LECT(NGFACE)
266 : ELSE
267 : NLFRI=0
268 : ENDIF
269 : IF (NLFRI.GT.0) THEN
270 : RHOF=MPRIMP.VPOCHA(NLFRI,1)
271 : ELSE
272 : RHOG = MPROC.VPOCHA(NLCGAU,1)
273 : RHOD = MPROC.VPOCHA(NLCDRO,1)
274 : RHOF = UMALPH * RHOG + ALPHA * RHOD
275 : ENDIF
276 : DKG = MPDKC.VPOCHA(NLCGAU,1)
277 : DKD = MPDKC.VPOCHA(NLCDRO,1)
278 : DKF = UMALPH * DKG + ALPHA * DKD
279 : NPP=NPTEL-1
280 : * IPD=1 : point à gauche du point NGFACE
281 : * IPD=2 : point à droite du point NGFACE
282 : DO 122 IPD=1,NPD
283 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
284 : IF (.NOT.LMUR) THEN
285 : CTGRY.POIDU2(IPD,IEL2)=NPDUAL
286 : ELSE
287 : CTGRY.POIDU1(IPD,IEL1)=NPDUAL
288 : ENDIF
289 : NLCEND=KRCENT.LECT(NPDUAL)
290 : IF (NLCEND.EQ.0) THEN
291 : WRITE(IOIMP,*) 'Erreur grave n°1'
292 : GOTO 9999
293 : ENDIF
294 : DO 124 IPP=1,NPP
295 : NPPRIM=MCOGRY.NUM(IPP+1,IELEM)
296 : LCTRB2=.TRUE.
297 : IF (LCLIMY) THEN
298 : NLCLY=KRYIMP.LECT(NPPRIM)
299 : IF (NLCLY.NE.0) THEN
300 : LCTRB2=.FALSE.
301 : ENDIF
302 : ENDIF
303 : IF (.NOT.LCTRB2) THEN
304 : * Lorsque une contribution est nulle, on fixe artificiellement le
305 : * point primal égal au point dual.
306 : IF (.NOT.LMUR) THEN
307 : CTGRY.POIPR2(IPP,IEL2)=NPDUAL
308 : RYKRHO.VALMA2(IPD,IPP,IEL2)=0.D0
309 : RYKRYK.VALMA2(IPD,IPP,IEL2)=0.D0
310 : ELSE
311 : CTGRY.POIPR1(IPP,IEL1)=NPDUAL
312 : RYKRHO.VALMA1(IPD,IPP,IEL1)=0.D0
313 : RYKRYK.VALMA1(IPD,IPP,IEL1)=0.D0
314 : ENDIF
315 : ELSE
316 : * Les contributions valent :
317 : * (d Res_{\rho Yk})_d / (d var)_p =
318 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \rho_f * Dk_f
319 : * * [ ((n_x * \beta_x) +(n_y * \beta_y) +(n_z * \beta_z)) *
320 : * ((dYk)_p / (d var)_p)]
321 : * avec :
322 : * (dYk)_p / (d \rho)_p = - Yk_p / \rho_p
323 : * (dYk)_p / (d \rho Yk)_p = 1 / \rho_p
324 : * \beta_x : coefficients pour le calcul de dYk/dx
325 : * \beta_y : coefficients pour le calcul de dYk/dy
326 : * \beta_z : coefficients pour le calcul de dYk/dz
327 : *
328 : NLCENP=KRCENT.LECT(NPPRIM)
329 : IF (NLCENP.EQ.0) THEN
330 : WRITE(IOIMP,*) 'Erreur grave n°2'
331 : GOTO 9999
332 : ENDIF
333 : * normale sortante pour IPD=1, rentrante pour IPD=2
334 : SIGNOR=(-1.D0)**(IPD+1)
335 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
336 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
337 : CNX =MPNORM.VPOCHA(NLFACE,1)
338 : CNY =MPNORM.VPOCHA(NLFACE,2)
339 : CNZ =MPNORM.VPOCHA(NLFACE,3)
340 : BETAX =KDYKDX.VELCHE(IPP+1,IELEM)
341 : BETAY =KDYKDY.VELCHE(IPP+1,IELEM)
342 : BETAZ =KDYKDZ.VELCHE(IPP+1,IELEM)
343 : RHOP =MPROC.VPOCHA(NLCENP,1)
344 : YKP =MPYKC.VPOCHA(NLCENP,1)
345 : * YKP =MPYKC.VPOCHA(NLCEND,1)
346 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*RHOF*DKF
347 : $ *((CNX*BETAX)+(CNY*BETAY)+(CNZ*BETAZ))
348 : DYKDRO=-YKP/RHOP
349 : DYKDRY=1.D0/RHOP
350 : IF (.NOT.LMUR) THEN
351 : CTGRY.POIPR2(IPP,IEL2)=NPPRIM
352 : RYKRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRO
353 : RYKRYK.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRY
354 : ELSE
355 : CTGRY.POIPR1(IPP,IEL1)=NPPRIM
356 : RYKRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRO
357 : RYKRYK.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRY
358 : ENDIF
359 : ENDIF
360 : 124 CONTINUE
361 : 122 CONTINUE
362 : IF (.NOT.LMUR) THEN
363 : IEL2=IEL2+1
364 : ELSE
365 : IEL1=IEL1+1
366 : ENDIF
367 : 12 CONTINUE
368 : NPP1=NPTEL-1
369 : NPP2=NPTEL-1
370 : NEL1=IEL1-1
371 : NEL2=IEL2-1
372 : SEGADJ RYKRHO
373 : SEGADJ RYKRYK
374 : SEGADJ CTGRY
375 : CTGRY.LMATSI(**)=RYKRHO
376 : CTGRY.LMATSI(**)=RYKRYK
377 : * On accumule les matrices résultantes dans IJACO
377 : CALL AJMTK(CTGRY,IJACO,IMPR,IRET)
378 : IF (IRET.NE.0) GOTO 9999
379 : SEGSUP RYKRHO
380 : SEGSUP RYKRYK
381 : SEGSUP CTGRY
382 : *
383 : SEGDES MCOGRY
384 : SEGDES KDYKDZ
385 : SEGDES KDYKDY
386 : SEGDES KDYKDX
387 : 1 CONTINUE
388 : SEGDES ICOGRY
389 : SEGDES MPYKC
390 : SEGDES MPROC
391 : SEGDES MPDKC
392 : SEGDES MPVOLU
393 : SEGDES MPNORM
394 : SEGDES MPSURF
395 : SEGDES MELEFL
396 : SEGDES KRFACE
397 : SEGDES KRCENT
398 : SEGDES NOMINC
399 : IF (LCLIMY) THEN
400 : SEGDES KRYIMP
401 : ENDIF
402 : IF (LCLIMR) THEN
403 : SEGDES KRRIMP
404 : SEGDES MPRIMP
405 : ENDIF
406 : ENDDO
407 : SEGDES PROPH2
408 : SEGDES PROPHY
409 : *
410 : * Normal termination
411 : *
412 : IRET=0
413 : RETURN
414 : *
415 : * Format handling
416 : *
417 : *
418 : * Error handling
419 : *
420 : 9999 CONTINUE
421 : IRET=1
422 : WRITE(IOIMP,*) 'An error was detected in subroutine zlap2b'
423 : RETURN
424 : *
425 : * End of subroutine ZLAP2B
426 : *
427 : END
428 :
429 :
430 :
431 :
432 :
433 :
434 :
435 :
© Cast3M 2003 - All rights reserved.
Disclaimer