1 : C ZLAP1B SOURCE MAGN 09/07/30 21:23:37 6458
2 : SUBROUTINE ZLAP1B(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 : ZLAP1B
12 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
13 : C VF 2D (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 : ZLAP1A : Calcul de la matrice jacobienne du
31 : C résidu du laplacien VF 2D.
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 2D.
55 : C SORTIES : -
56 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
57 : C***********************************************************************
58 : C VERSION : v1, 22/02/2002, version initiale
59 : C HISTORIQUE : v1, 22/02/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
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,CNX,CNY
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,XG,YG,XFMXG,YFMYG,DRG
148 : & ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
149 : C
150 : *
151 : * Executable statements
152 : *
153 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans zlap1b.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 : SEGDES JCOGRY
189 : SEGACT KDYKDX
190 : SEGACT KDYKDY
191 : SEGACT MCOGRY
192 : NELEM=MCOGRY.NUM(/2)
193 : NPTEL=MCOGRY.NUM(/1)
194 : NPP1=NPTEL-1
195 : NPP2=NPTEL-1
196 : NEL1=NELEM
197 : NEL2=NELEM
198 : IEL1=1
199 : IEL2=1
200 : SEGINI RYKRHO
201 : SEGINI RYKRYK
202 : SEGINI CTGRY
203 : RYKRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
204 : RYKRHO.NOMPRI(5:8)=' '
205 : RYKRHO.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
206 : RYKRHO.NOMDUA(5:8)=' '
207 : RYKRYK.NOMPRI(1:4)=NOMINC.MOTS(IDIM+2+IESP)
208 : RYKRYK.NOMPRI(5:8)=' '
209 : RYKRYK.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
210 : RYKRYK.NOMDUA(5:8)=' '
211 : DO 12 IELEM=1,NELEM
212 : * Le premier point du support de ICOGRY est un point FACE
213 : NGFACE=MCOGRY.NUM(1,IELEM)
214 : NLFACE=KRFACE.LECT(NGFACE)
215 : IF (NLFACE.EQ.0) THEN
216 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
217 : GOTO 9999
218 : ENDIF
219 : * On calcule la contribution à la matrice jacobienne IJACO de la face
220 : * NGFAC (points duaux : centres à gauche et à droite de la face)
221 : * (points primaux : une partie (bicoz conditions aux limites)
222 : * de ceux du stencil pour le calcul du gradient
223 : * à la face, ils doivent être des points centres)
224 : * Si le flux de chaleur sur la face est imposé par les conditions
225 : * aux limites, la contribution de la face à IJACO est nulle.
226 : NGCGAU=MELEFL.NUM(1,NLFACE)
227 : NGCDRO=MELEFL.NUM(3,NLFACE)
228 : NLCGAU=KRCENT.LECT(NGCGAU)
229 : NLCDRO=KRCENT.LECT(NGCDRO)
230 : LMUR=(NGCGAU.EQ.NGCDRO)
231 : * On distingue le cas où la face est un bord du maillage (mur)
232 : * du cas où la face est interne au maillage
233 : IF (.NOT.LMUR) THEN
234 : NPD=2
235 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
236 : XF = MCOORD.XCOOR(ICOORX)
237 : YF = MCOORD.XCOOR(ICOORX+1)
238 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
239 : XG = MCOORD.XCOOR(ICOORX)
240 : YG = MCOORD.XCOOR(ICOORX+1)
241 : XFMXG = XF - XG
242 : YFMYG = YF - YG
243 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
244 : ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
245 : XD = MCOORD.XCOOR(ICOORX)
246 : YD = MCOORD.XCOOR(ICOORX+1)
247 : XFMXD = XF - XD
248 : YFMYD = YF - YD
249 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
250 : ALPHA=DRG/(DRG+DRD)
251 : UMALPH= 1.0D0 - ALPHA
252 : ELSE
253 : NPD=1
254 : ALPHA=0.0D0
255 : UMALPH=1.0D0
256 : ENDIF
257 : IF (LCLIMR) THEN
258 : NLFRI=KRRIMP.LECT(NGFACE)
259 : ELSE
260 : NLFRI=0
261 : ENDIF
262 : IF (NLFRI.GT.0) THEN
263 : RHOF=MPRIMP.VPOCHA(NLFRI,1)
264 : ELSE
265 : RHOG = MPROC.VPOCHA(NLCGAU,1)
266 : RHOD = MPROC.VPOCHA(NLCDRO,1)
267 : RHOF = UMALPH * RHOG + ALPHA * RHOD
268 : ENDIF
269 : DKG = MPDKC.VPOCHA(NLCGAU,1)
270 : DKD = MPDKC.VPOCHA(NLCDRO,1)
271 : DKF = UMALPH * DKG + ALPHA * DKD
272 : NPP=NPTEL-1
273 : * IPD=1 : point à gauche du point NGFACE
274 : * IPD=2 : point à droite du point NGFACE
275 : DO 122 IPD=1,NPD
276 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
277 : IF (.NOT.LMUR) THEN
278 : CTGRY.POIDU2(IPD,IEL2)=NPDUAL
279 : ELSE
280 : CTGRY.POIDU1(IPD,IEL1)=NPDUAL
281 : ENDIF
282 : NLCEND=KRCENT.LECT(NPDUAL)
283 : IF (NLCEND.EQ.0) THEN
284 : WRITE(IOIMP,*) 'Erreur grave n°1'
285 : GOTO 9999
286 : ENDIF
287 : DO 124 IPP=1,NPP
288 : NPPRIM=MCOGRY.NUM(IPP+1,IELEM)
289 : LCTRB2=.TRUE.
290 : IF (LCLIMY) THEN
291 : NLCLY=KRYIMP.LECT(NPPRIM)
292 : IF (NLCLY.NE.0) THEN
293 : LCTRB2=.FALSE.
294 : ENDIF
295 : ENDIF
296 : IF (.NOT.LCTRB2) THEN
297 : * Lorsque une contribution est nulle, on fixe artificiellement le
298 : * point primal égal au point dual.
299 : IF (.NOT.LMUR) THEN
300 : CTGRY.POIPR2(IPP,IEL2)=NPDUAL
301 : RYKRHO.VALMA2(IPD,IPP,IEL2)=0.D0
302 : RYKRYK.VALMA2(IPD,IPP,IEL2)=0.D0
303 : ELSE
304 : CTGRY.POIPR1(IPP,IEL1)=NPDUAL
305 : RYKRHO.VALMA1(IPD,IPP,IEL1)=0.D0
306 : RYKRYK.VALMA1(IPD,IPP,IEL1)=0.D0
307 : ENDIF
308 : ELSE
309 : * Les contributions valent :
310 : * (d Res_{\rho Yk})_d / (d var)_p =
311 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \rho_f * Dk_f
312 : * * [ ((n_x * \beta_x) +(n_y * \beta_y)) *
313 : * ((dYk)_p / (d var)_p)]
314 : * avec :
315 : * (dYk)_p / (d \rho)_p = - Yk_p / \rho_p
316 : * (dYk)_p / (d \rho Yk)_p = 1 / \rho_p
317 : * \beta_x : coefficients pour le calcul de dYk/dx
318 : * \beta_y : coefficients pour le calcul de dYk/dy
319 : *
320 : NLCENP=KRCENT.LECT(NPPRIM)
321 : IF (NLCENP.EQ.0) THEN
322 : WRITE(IOIMP,*) 'Erreur grave n°2'
323 : GOTO 9999
324 : ENDIF
325 : * normale sortante pour IPD=1, rentrante pour IPD=2
326 : SIGNOR=(-1.D0)**(IPD+1)
327 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
328 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
329 : CNX =MPNORM.VPOCHA(NLFACE,1)
330 : CNY =MPNORM.VPOCHA(NLFACE,2)
331 : BETAX =KDYKDX.VELCHE(IPP+1,IELEM)
332 : BETAY =KDYKDY.VELCHE(IPP+1,IELEM)
333 : RHOP =MPROC.VPOCHA(NLCENP,1)
334 : YKP =MPYKC.VPOCHA(NLCENP,1)
335 : * YKP =MPYKC.VPOCHA(NLCEND,1)
336 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*RHOF*DKF
337 : $ *((CNX*BETAX)+(CNY*BETAY))
338 : DYKDRO=-YKP/RHOP
339 : DYKDRY=1.D0/RHOP
340 : IF (.NOT.LMUR) THEN
341 : CTGRY.POIPR2(IPP,IEL2)=NPPRIM
342 : RYKRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRO
343 : RYKRYK.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRY
344 : ELSE
345 : CTGRY.POIPR1(IPP,IEL1)=NPPRIM
346 : RYKRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRO
347 : RYKRYK.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRY
348 : ENDIF
349 : ENDIF
350 : 124 CONTINUE
351 : 122 CONTINUE
352 : IF (.NOT.LMUR) THEN
353 : IEL2=IEL2+1
354 : ELSE
355 : IEL1=IEL1+1
356 : ENDIF
357 : 12 CONTINUE
358 : NPP1=NPTEL-1
359 : NPP2=NPTEL-1
360 : NEL1=IEL1-1
361 : NEL2=IEL2-1
362 : SEGADJ RYKRHO
363 : SEGADJ RYKRYK
364 : SEGADJ CTGRY
365 : CTGRY.LMATSI(**)=RYKRHO
366 : CTGRY.LMATSI(**)=RYKRYK
367 : * On accumule les matrices résultantes dans IJACO
367 : CALL AJMTK(CTGRY,IJACO,IMPR,IRET)
368 : IF (IRET.NE.0) GOTO 9999
369 : SEGSUP RYKRHO
370 : SEGSUP RYKRYK
371 : SEGSUP CTGRY
372 : *
373 : SEGDES MCOGRY
374 : SEGDES KDYKDY
375 : SEGDES KDYKDX
376 : 1 CONTINUE
377 : SEGDES ICOGRY
378 : SEGDES MPYKC
379 : SEGDES MPROC
380 : SEGDES MPDKC
381 : SEGDES MPVOLU
382 : SEGDES MPNORM
383 : SEGDES MPSURF
384 : SEGDES MELEFL
385 : SEGDES KRFACE
386 : SEGDES KRCENT
387 : SEGDES NOMINC
388 : IF (LCLIMY) THEN
389 : SEGDES KRYIMP
390 : ENDIF
391 : IF (LCLIMR) THEN
392 : SEGDES KRRIMP
393 : SEGDES MPRIMP
394 : ENDIF
395 : ENDDO
396 : SEGDES PROPH2
397 : SEGDES PROPHY
398 : *
399 : * Normal termination
400 : *
401 : IRET=0
402 : RETURN
403 : *
404 : * Format handling
405 : *
406 : *
407 : * Error handling
408 : *
409 : 9999 CONTINUE
410 : IRET=1
411 : WRITE(IOIMP,*) 'An error was detected in subroutine zlap1b'
412 : RETURN
413 : *
414 : * End of subroutine ZLAP1B
415 : *
416 : END
417 :
418 :
419 :
420 :
421 :
422 :
423 :
424 :
© Cast3M 2003 - All rights reserved.
Disclaimer