1 : C KONMS3 SOURCE MAGN 09/07/30 21:18:27 6458
2 : SUBROUTINE KONMS3(NSP,ILIINC,IRN,IUN,IPN,IYN,MLRECP,
3 : & MLRECV,INORM,ICHPVO,ICHPSU,MELEMC,MELEFE,MELLIM,IMAT)
4 : C
5 : C************************************************************************
6 : C
7 : C PROJET : CASTEM 2000
8 : C
9 : C NOM : KONMS3 ("convection for multispecies")
10 : C
11 : C DESCRIPTION : Voir KON19 (appele par KON19)
12 : C Calcul du jacobien du résidu pour la méthode
13 : C AUSM+
14 : C
14 : C Cas trois dimensions, gaz "calorically perfect"
15 : C MULTISPECIES!!!!!!!!
16 : C
17 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
18 : C
19 : C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
20 : C
21 : C************************************************************************
22 : C
23 : C
24 : C APPELES (Outils
25 : C CASTEM) : KRIPAD, LICHT, ERREUR
26 : C
27 : C APPELES (Calcul) : CMS3D, CWMS3D
28 : C
29 : C************************************************************************
30 : C
31 : C ENTREES
32 : C
33 : C NSP : number of species (total) ;
34 : C
35 : C ILINC : liste des inconnues (pointeur d'un objet de type LISTMOTS)
36 : C
37 : C 1) Pointeurs des CHPOINT
38 : C
39 : C IRN : CHPOINT CENTRE contenant la masse volumique ;
40 : C
41 : C IUN : CHPOINT CENTRE contenant la vitesse ;
42 : C
43 : C IPN : CHPOINT CENTRE contenant la pression ;
44 : C
45 : C IYN : CHPOINT CENTRE contenant les fractions massiques ;
46 : C
47 : C INORM : CHPOINT FACE contenant les normales aux faces ;
48 : C
49 : C ICHPVO : CHPOINT VOLUME contenant le volume
50 : C
51 : C ICHPSU : CHPOINT FACE contenant la surface des faces
52 : C
53 : C
54 : C 2) Pointeurs des LIST REELS
55 : C
56 : C MLRECP : LIST REELS contenant les CP's des gases differents ;
57 : C
58 : C MLRECV : LIST REELS contenant les CV's des gases differents ;
59 : C
60 : C 3) Pointeurs de MELEME de la table DOMAINE
61 : C
62 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
63 : C
64 : C MELEFE : MELEME 'FACEL' du connectivité Faces -> Elts
65 : C
66 : C MELLIM : MELEME SPG des conditions aux bords
67 : C
68 : C SORTIES
69 : C
70 : C IMAT : pointeur de la MATRIK du jacobien du residu
71 : C
72 : C************************************************************************
73 : C
74 : C HISTORIQUE (Anomalies et modifications éventuelles)
75 : C
76 : C HISTORIQUE :
77 : C
78 : C************************************************************************
79 : C
80 : C
81 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0
82 : C GAMMA \in (1,3)
83 : C Si non il faut le faire!!!
84 : C
85 : C************************************************************************
86 : C
87 : C
88 : C**** Variables de COOPTIO
89 : C
90 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
91 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
92 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
93 : C & ,IECHO, IIMPI, IOSPI
94 : C & ,IDIM
95 : C & ,MCOORD
96 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
97 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
98 : C & ,NORINC,NORVAL,NORIND,NORVAD
99 : C & ,NUCROU, IPSAUV
100 : C---------------------------------------------------
101 : IMPLICIT INTEGER(I-N)
102 : INTEGER NSP,II,JJ,KV,I,J,IFICLE,IPREFI,JLL,JRR
103 : INTEGER ILIINC, IRN,IUN,IPN,IYN,INORM,ICHPVO,ICHPSU
104 : & , IMAT, IGEOMC, IGEOMF
105 : & , NFAC, NBSOUS, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID
106 : & , NKMT, NBME, NBEL, MP, NP
107 : & , IFAC, NGCF, NLCF, NGCG, NGCD, NLCG, NLCD, NLFL
108 : REAL*8 ROG, PG, UXG, UYG, UZG, VOLG
109 : & , ROD, PD, UXD, UYD, UZD, VOLD
110 : & , SURF, FUNCEL
111 : REAL*8 WVEC_L(5), WVEC_R(5), NVECT(3), TVECT1(3), TVECT2(3)
112 : REAL*8 C11,C12,C13,C21,C22,C23,C31,C32,C33,DET
113 : REAL*8 ZC11,ZC12,ZC13,ZC21,ZC22,ZC23,ZC31,ZC32,ZC33
114 : CHARACTER*8 TYPE
115 : C
116 : C**** LES INCLUDES
117 : C
118 : -INC CCOPTIO
119 : -INC SMCHPOI
120 : -INC SMMATRIK
121 : -INC SMELEME
122 : -INC SMLMOTS
123 : -INC SMLENTI
124 : POINTEUR MPRN.MPOVAL, MPUN.MPOVAL, MPPN.MPOVAL, MPYN.MPOVAL,
125 : & MPNORM.MPOVAL, MPVOLU.MPOVAL, MPOVSU.MPOVAL
126 : POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME,
127 : & MELEDU.MELEME, MELLIM.MELEME
128 : POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLELIM.MLENTI
129 : POINTEUR CELL.IZAFM
130 : POINTEUR MLMINC.MLMOTS
131 : C----------------------------------------------------
132 : -INC SMLREEL
133 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
134 : C-------------------------------------------------------
135 : C********* Les Jacobians ******************************
136 : C-------------------------------------------------------
137 : SEGMENT JACEL
138 : REAL*8 JAC(4+NSP,4+NSP)
139 : ENDSEGMENT
140 : POINTEUR JTL.JACEL, JTR.JACEL, JTT.JACEL
141 : C--------------------------------------------------------
142 : C KRIPAD pour la correspondance global/local des centres
143 : C--------------------------------------------------------
143 : CALL KRIPAD(MELLIM,MLELIM)
143 : CALL KRIPAD(MELEMC,MLENTC)
144 : C------------------
145 : SEGACT MELEMC
146 : C------------------
147 : SEGACT MELEFE
148 : C------------------
148 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
148 : CALL LICHT(INORM,MPNORM,TYPE,IGEOMF)
148 : CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC)
149 : C----------------------------------------------
150 : MELEMF = IGEOMF
150 : CALL KRIPAD(MELEMF,MLENTF)
151 : SEGACT MELEMF
152 : C-----------------------------------------------
152 : CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
152 : CALL LICHT(IPN,MPPN,TYPE,IGEOMC)
152 : CALL LICHT(IUN,MPUN,TYPE,IGEOMC)
152 : CALL LICHT(IYN,MPYN,TYPE,IGEOMC)
153 : C-----------------------------------------------
154 : NFAC = MELEFE.NUM(/2)
155 : C-----------------------------------------------
156 : C**** Maillage des inconnues primales
157 : C-----------------------------------------------
158 : NBSOUS = 0
159 : NBREF = 0
160 : NBELEM = NFAC
161 : NBNN = 2
162 : C-----------------------------------------------
163 : SEGINI MELEDU
164 : C----------------------
165 : C** MELEDU = 'SEG2' **
166 : C----------------------
167 : MELEDU.ITYPEL = 2
168 : C--------------
169 : NRIGE = 7
170 : NMATRI = 1
171 : NKID = 9
172 : NKMT = 7
173 : C--------------
174 : SEGINI MATRIK
175 : IMAT = MATRIK
176 : MATRIK.IRIGEL(1,1) = MELEDU
177 : MATRIK.IRIGEL(2,1) = MELEDU
178 : C-------------------------
179 : C Matrice non symetrique
180 : C-------------------------
181 : MATRIK.IRIGEL(7,1) = 2
182 : C-------------------------
183 : NBME = (4+NSP)*(4+NSP)
184 : NBSOUS = 1
185 : SEGINI IMATRI
186 : MLMINC = ILIINC
187 : SEGACT MLMINC
188 : MATRIK.IRIGEL(4,1) = IMATRI
189 : C-----------------------------------------------
190 : DO 1 J=1,(NSP+4)
191 : KV=(J-1)*(4+NSP)
192 : IMATRI.LISPRI(KV+1) = MLMINC.MOTS(1)
193 : IMATRI.LISPRI(KV+2) = MLMINC.MOTS(2)
194 : IMATRI.LISPRI(KV+3) = MLMINC.MOTS(3)
195 : IMATRI.LISPRI(KV+4) = MLMINC.MOTS(4)
196 : IMATRI.LISPRI(KV+5) = MLMINC.MOTS(5)
197 : DO 2 I=1,(NSP-1)
198 : IMATRI.LISPRI(KV+5+I) = MLMINC.MOTS(5+I)
199 : 2 CONTINUE
200 : 1 CONTINUE
201 : C-----------------------------------------------
202 : DO 3 J=1,(NSP+4)
203 : KV=(J-1)*(4+NSP)
204 : IMATRI.LISDUA(KV+1) = MLMINC.MOTS(j)
205 : IMATRI.LISDUA(KV+2) = MLMINC.MOTS(j)
206 : IMATRI.LISDUA(KV+3) = MLMINC.MOTS(j)
207 : IMATRI.LISDUA(KV+4) = MLMINC.MOTS(j)
208 : IMATRI.LISDUA(KV+5) = MLMINC.MOTS(j)
209 : DO 4 I=1,(NSP-1)
210 : IMATRI.LISDUA(KV+5+I) = MLMINC.MOTS(j)
211 : 4 CONTINUE
212 : 3 CONTINUE
213 : C-----------------------------------------------
214 : C-----------------------------------------------
215 : NBEL = NBELEM
216 : NBSOUS = 1
217 : NP = 2
218 : MP = 2
219 : C-----------------------------------------------------------
220 : C-----------------------------------------------------------
221 : DO 5 I=1,NBME
222 : SEGINI CELL
223 : IMATRI.LIZAFM(1,I) = CELL
224 : 5 CONTINUE
225 : C---------------------------------
226 : C---------------------------------
227 : DO IFAC = 1, NFAC, 1
228 : NGCF = MELEFE.NUM(2,IFAC)
229 : NLCF = MLENTF.LECT(NGCF)
230 : IF(NLCF .NE. IFAC)THEN
231 : WRITE(IOIMP,*) 'Il ne faut pas jouer avec la table domaine'
231 : CALL ERREUR(5)
232 : GOTO 9999
233 : ENDIF
234 : NLFL = MLELIM.LECT(NGCF)
235 : NGCG = MELEFE.NUM(1,IFAC)
236 : NGCD = MELEFE.NUM(3,IFAC)
237 : IF(NLFL .NE. 0)THEN
238 : C--------------------------------------------------------
239 : C The point belongs on BC -> No contribution to jacobian!
240 : C--------------------------------------------------------
241 : MELEDU.NUM(1,IFAC) = NGCG
242 : MELEDU.NUM(2,IFAC) = NGCD
243 : ELSEIF(NGCG .NE. NGCD)THEN
244 : C-----------------------------------
245 : C********** Les MELEMEs
246 : C-----------------------------------
247 : MELEDU.NUM(1,IFAC) = NGCG
248 : MELEDU.NUM(2,IFAC) = NGCD
249 : C-----------------------------------
250 : C********** Les etats G et D
251 : C-----------------------------------
252 : NLCG = MLENTC.LECT(NGCG)
253 : NLCD = MLENTC.LECT(NGCD)
254 : C-----------------
255 : ROG = MPRN.VPOCHA(NLCG,1)
256 : PG = MPPN.VPOCHA(NLCG,1)
257 : UXG = MPUN.VPOCHA(NLCG,1)
258 : UYG = MPUN.VPOCHA(NLCG,2)
259 : UZG = MPUN.VPOCHA(NLCG,3)
260 : VOLG = MPVOLU.VPOCHA(NLCG,1)
261 : C-------------------------------------------------
262 : WVEC_L(1)=ROG
263 : WVEC_L(2)=UXG
264 : WVEC_L(3)=UYG
265 : WVEC_L(4)=UZG
266 : WVEC_L(5)=PG
267 : C-------------------------------------------------
268 : ROD = MPRN.VPOCHA(NLCD,1)
269 : PD = MPPN.VPOCHA(NLCD,1)
270 : UXD = MPUN.VPOCHA(NLCD,1)
271 : UYD = MPUN.VPOCHA(NLCD,2)
272 : UZD = MPUN.VPOCHA(NLCD,3)
273 : VOLD = MPVOLU.VPOCHA(NLCD,1)
274 : C------------------------------------------------
275 : WVEC_R(1)=ROD
276 : WVEC_R(2)=UXD
277 : WVEC_R(3)=UYD
278 : WVEC_R(4)=UZD
279 : WVEC_R(5)=PD
280 : C------------------------------------------------
281 : C La normale G->D
282 : C La tangente
283 : C------------------------------------------------
284 : SURF = MPOVSU.VPOCHA(NLCF,1)
285 : NVECT(1) = MPNORM.VPOCHA(NLCF,7)
286 : NVECT(2) = MPNORM.VPOCHA(NLCF,8)
287 : NVECT(3) = MPNORM.VPOCHA(NLCF,9)
288 : c-----------------------------------------------
289 : TVECT1(1) = MPNORM.VPOCHA(NLCF,1)
290 : TVECT1(2) = MPNORM.VPOCHA(NLCF,2)
291 : TVECT1(3) = MPNORM.VPOCHA(NLCF,3)
292 : c----------------------------------------------
293 : TVECT2(1) = MPNORM.VPOCHA(NLCF,4)
294 : TVECT2(2) = MPNORM.VPOCHA(NLCF,5)
295 : TVECT2(3) = MPNORM.VPOCHA(NLCF,6)
296 : C-------------------------------------------------
296 : CALL cms3d(NSP,JLL,JRR,WVEC_L,WVEC_R,NVECT,TVECT1,TVECT2,
297 : & MPYN,MLRECP,MLRECV,NLCG,NLCD)
298 : C-----------------------------------------------------
299 : C********** AB.AM(IFAC,IPRIM,IDUAL)
300 : C A = nom de l'inconnu duale (Ro,rUX,rUY,RET)
301 : C B = nom de l'inconnu primale (Ro,rUX,rUY,RET)
302 : C IPRIM = 1, 2 -> G, D
303 : C IDUAL = 1, 2 -> G, D
304 : C i.e.
305 : C A_IDUAL = AB.AM(IFAC,IPRIM,IDUAL) * B_IPRIM + ...
306 : C
307 : C
308 : C********** Dual RN
309 : C----------------------------------------------------
310 : JTL=JLL
311 : JTR=JRR
312 : SEGACT JTL
313 : SEGACT JTR
314 : DO 10 II = 1,(4+NSP)
315 : DO 15 JJ = 1,(4+NSP)
316 : KV = (II-1)*(4+NSP)
317 : C----------------------------------
318 : CELL = IMATRI.LIZAFM(1,KV+JJ)
319 : FUNCEL = SURF * JTL.JAC(II,JJ)
320 : CELL.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
321 : CELL.AM(IFAC,1,2) = FUNCEL / VOLD
322 : 15 CONTINUE
323 : 10 CONTINUE
324 : C-----------------------------------------------------
325 : DO 20 II = 1,(4+NSP)
326 : DO 25 JJ = 1,(4+NSP)
327 : KV = (II-1)*(4+NSP)
328 : CELL = IMATRI.LIZAFM(1,KV+JJ)
329 : FUNCEL = SURF * JTR.JAC(II,JJ)
330 : CELL.AM(IFAC,2,2) = FUNCEL / VOLD
331 : CELL.AM(IFAC,2,1) = -FUNCEL / VOLG
332 : 25 CONTINUE
333 : 20 CONTINUE
334 : SEGDES JTR
335 : SEGDES JTL
336 : C-----------------------------------------------------
337 : ELSE
338 : C-----------------------------------------------------
339 : C************* Murs (NGCG = NGCD) ******************
340 : C-----------------------------------------------------
341 : MELEDU.NUM(1,IFAC) = NGCG
342 : MELEDU.NUM(2,IFAC) = NGCD
343 : NLCG = MLENTC.LECT(NGCG)
344 : C--------------------------------------
345 : ROG = MPRN.VPOCHA(NLCG,1)
346 : PG = MPPN.VPOCHA(NLCG,1)
347 : UXG = MPUN.VPOCHA(NLCG,1)
348 : UYG = MPUN.VPOCHA(NLCG,2)
349 : UZG = MPUN.VPOCHA(NLCG,3)
350 : VOLG = MPVOLU.VPOCHA(NLCG,1)
351 : C-------------------------------------------
352 : WVEC_L(1)=ROG
353 : WVEC_L(2)=UXG
354 : WVEC_L(3)=UYG
355 : WVEC_L(4)=UZG
356 : WVEC_L(5)=PG
357 : C-------------------------------------------------
358 : SURF = MPOVSU.VPOCHA(NLCF,1)
359 : NVECT(1) = MPNORM.VPOCHA(NLCF,7)
360 : NVECT(2) = MPNORM.VPOCHA(NLCF,8)
361 : NVECT(3) = MPNORM.VPOCHA(NLCF,9)
362 : c--------------------------------------------
363 : TVECT1(1) = MPNORM.VPOCHA(NLCF,1)
364 : TVECT1(2) = MPNORM.VPOCHA(NLCF,2)
365 : TVECT1(3) = MPNORM.VPOCHA(NLCF,3)
366 : c----------------------------------------------
367 : TVECT2(1) = MPNORM.VPOCHA(NLCF,4)
368 : TVECT2(2) = MPNORM.VPOCHA(NLCF,5)
369 : TVECT2(3) = MPNORM.VPOCHA(NLCF,6)
370 : C------- COEFFICIENTS ----------------------------
371 : C11=TVECT1(2)*TVECT2(3)-TVECT1(3)*TVECT2(2)
372 : C12=NVECT(2)*TVECT2(3)-TVECT2(2)*NVECT(3)
373 : C13=NVECT(2)*TVECT1(3)-TVECT1(2)*NVECT(3)
374 : C---------------------------------
375 : C21=TVECT1(1)*TVECT2(3)-TVECT1(3)*TVECT2(1)
376 : C22=NVECT(1)*TVECT2(3)-TVECT2(1)*NVECT(3)
377 : C23=NVECT(1)*TVECT1(3)-TVECT1(1)*NVECT(3)
378 : C---------------------------------
379 : C31=TVECT1(1)*TVECT2(2)-TVECT1(2)*TVECT2(1)
380 : C32=NVECT(1)*TVECT2(2)-TVECT2(1)*NVECT(2)
381 : C33=NVECT(1)*TVECT1(2)-TVECT1(1)*NVECT(2)
382 : DET=NVECT(1)*C11-NVECT(2)*C21+NVECT(3)*C31
383 : C---------------------------------
384 : ZC11=-NVECT(1)*C11-TVECT1(1)*C12+TVECT2(1)*C13
385 : ZC12=-NVECT(2)*C11-TVECT1(2)*C12+TVECT2(2)*C13
386 : ZC13=-NVECT(3)*C11-TVECT1(3)*C12+TVECT2(3)*C13
387 : C---------------------------------
388 : ZC21=NVECT(1)*C21+TVECT1(1)*C22-TVECT2(1)*C23
389 : ZC22=NVECT(2)*C21+TVECT1(2)*C22-TVECT2(2)*C23
390 : ZC23=NVECT(3)*C21+TVECT1(3)*C22-TVECT2(3)*C23
391 : C---------------------------------
392 : ZC31=-NVECT(1)*C31-TVECT1(1)*C32+TVECT2(1)*C33
393 : ZC32=-NVECT(2)*C31-TVECT1(2)*C32+TVECT2(2)*C33
394 : ZC33=-NVECT(3)*C31-TVECT1(3)*C32+TVECT2(3)*C33
395 : C-------------------------------------------------
396 : ROD = ROG
397 : PD = PG
398 : UXD = (ZC11*UXG+ZC12*UYG+ZC13*UZG)/DET
399 : UYD = (ZC21*UXG+ZC22*UYG+ZC23*UZG)/DET
400 : UZD = (ZC31*UXG+ZC32*UYG+ZC33*UZG)/DET
401 : VOLD = VOLG
402 : C------------------------------------------------
403 : WVEC_R(1)=ROD
404 : WVEC_R(2)=UXD
405 : WVEC_R(3)=UYD
406 : WVEC_R(4)=UZD
407 : WVEC_R(5)=PD
408 : C-------------------------------------------
409 : C********** La normale sortante
410 : C-------------------------------------------
410 : CALL cwms3d(NSP,JLL,WVEC_L,WVEC_R,NVECT,
411 : & TVECT1,TVECT2,MPYN,MLRECP,MLRECV,NLCG)
412 : C----------------------------------------------------
413 : JTL=JLL
414 : SEGACT JTL
415 : DO 70 II = 1,(4+NSP)
416 : DO 75 JJ = 1,(4+NSP)
417 : KV = (II-1)*(4+NSP)
418 : C---------------------------------
419 : CELL = IMATRI.LIZAFM(1,KV+JJ)
420 : FUNCEL = SURF * JTL.JAC(II,JJ)
421 : CELL.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
422 : CELL.AM(IFAC,1,2) = 0.0D0
423 : 75 CONTINUE
424 : 70 CONTINUE
425 : SEGDES JTL
426 : C--------------------------------------------
427 : DO 40 II = 1,(4+NSP)
428 : DO 45 JJ = 1,(4+NSP)
429 : KV = (II-1)*(4+NSP)
430 : CELL = IMATRI.LIZAFM(1,KV+JJ)
431 : CELL.AM(IFAC,2,2) = 0.0D0
432 : CELL.AM(IFAC,2,1) = 0.0D0
433 : 45 CONTINUE
434 : 40 CONTINUE
435 : C--------------------------------------------
436 : ENDIF
437 : ENDDO
438 : C
439 : SEGDES MELEMC
440 : SEGDES MELEFE
441 : SEGDES MELEMF
442 : C
443 : SEGDES MPOVSU
444 : SEGDES MPVOLU
445 : SEGDES MPNORM
446 : C
447 : SEGDES MPRN
448 : SEGDES MPPN
449 : SEGDES MPUN
450 : SEGDES MPYN
451 : C
452 : SEGDES MELEDU
453 : SEGDES MATRIK
454 : DO 80 II=1,NBME
455 : CELL = IMATRI.LIZAFM(1,II)
456 : SEGDES CELL
457 : 80 CONTINUE
458 : SEGDES IMATRI
459 : C
460 : SEGSUP MLENTC
461 : SEGSUP MLENTF
462 : SEGSUP MLELIM
463 : SEGDES MLMINC
464 : IF(MELLIM .NE.0) SEGDES MELLIM
465 :
466 : 9999 CONTINUE
467 : RETURN
468 : END
469 :
470 :
471 :
472 :
473 :
474 :
475 :
476 :
477 :
478 :
479 :
480 :
481 :
© Cast3M 2003 - All rights reserved.
Disclaimer