1 : C KFM1 SOURCE CHAT 06/08/24 21:46:35 5529
2 : SUBROUTINE KFM1
3 : C************************************************************************
4 : C
5 : C PROJET : CASTEM 2000
6 : C
7 : C NOM : KFM1
8 : C
9 : C DESCRIPTION : Subroutine appellée par KONV1
10 : C
11 : C Modelisation 2D/3D des equations d'Euler
12 : C
13 : C Free matrix method: calcul d'un stationaire
14 : C avec un jacobi par point, un jacobi par ligne,
15 : C un jacobi avec 2 sweeps
16 : C
17 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
18 : C
19 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF
20 : C
21 : C************************************************************************
22 : C
23 : C APPELES (Calcul) :
24 : C
25 : C************************************************************************
26 : C
27 : C*** SYNTAXE
28 : C
29 : C Discrétisation en VF "cell-centered" des équations d'Euler pour
30 : C un gaz parfait mono-constituent polytropique
31 : C Inconnues: densité, quantité de mouvement, énergie totale par
32 : C unité de volumes (variables conservatives UN)
33 : C
34 : C A(UN) '*' DUN_L = RES(UN) + B(UN) - B(UN_L)
35 : C
36 : C ou
37 : C
38 : C RES(UN) est calculé avec ('KONV' 'VF' 'PERFMONO')
39 : C UN_0 = UN
40 : C UN_L = UN + DUN_L
41 : C
42 : C Ici on calcule DU
43 : C
44 : C DUN = 'KONV' 'VF' 'PMON1FMM' 'JACOBI' LMOT1
45 : C MOD1 MCHPO0 MCHPO1 MCHPO2 MCHPO3 MCHPO4 FLOT0 FLOT1
46 : C ENT1 ('CLIM' LMOT2 MCHPO5 MCHPO7) MCHPO6 ;
47 : C
48 : C
49 : C LMOT1 : Il contient dans l'ordre suivant: le noms de la densité,
50 : C de la qdm, de l'énergie totale par unité de volume
51 : C
52 : C MOD1 : objet modele de type EULER
53 : C
54 : C MCHPO0 : CHPOINT contenant le residu explicite; composantes =
55 : C LMOT1, SPG = (DOMA MOD1 'CENTRE').
56 : C
57 : C MCHPO1 : CHPOINT contenant la masse volumique (en kg/m^3; une
58 : C composante, 'SCAL', SPG =(DOMA MOD1 'CENTRE')).
59 : C
60 : C MCHPO2 : CHPOINT contenant les débits (en kg/s/m^2; deux
61 : C composantes en 2D, 'UX ','UY ', trois composantes
62 : C en 3D, 'UX ','UY ', 'UZ ', SPG =(DOMA MOD1 'CENTRE')).
63 : C
64 : C MCHPO3 : CHPOINT contenant l'énergie totale par unité de volume
65 : C (en J/m^3; une composante, 'SCAL', SPG =(DOMA MOD1 'CENTRE')).
66 : C
67 : C MCHPO4 : CHPOINT contenant le "gamma" du gaz (une composante,
68 : C 'SCAL', SPG =(DOMA MOD1 'CENTRE')).
69 : C
70 : C LMOT2 : Il contient dans l'ordre suivant: le noms de la densité,
71 : C de la vitesse, de la pression.
72 : C
73 : C FLOT0 : pas de temps physique
74 : C
75 : C FLOT1 : facteur de securité (double de la CFL) pour le pas de temps
76 : C dual
77 : C
78 : C ENT1 : iterations dans le jacobi
79 : C
80 : C MCHPO5 : CHPOINT contenant les conditions aux limites
81 : C (composantes en LMOT2).
82 : C
83 : C MCHPO6 : CHPOINT contenant le coeff pour le calcul du rayon
84 : C spectral visc. (une composante,
85 : C 'SCAL', SPG =(DOMA MOD1 'CENTRE')).
86 : C
87 : C MCHPO7 : CHPOINT contenant les types de conditions aux limites
88 : C
89 : C************************************************************************
90 : C
91 : C HISTORIQUE (Anomalies et modifications éventuelles)
92 : C
93 : C HISTORIQUE : crée le 30/04/02
94 : C Janvier 2003: implementation de condition aux limites
95 : C
96 : C HISTORIQUE : rajouts par T. KLOCZKO entre 2003 et 2006
97 : C
98 : C Deux types de méthodes de relaxation disponibles
99 : C
100 : C - Point Jacobi ou PJ (kfm11.eso (2D) et kfm13.eso (3D))
101 : C - Symmetric Gauss-Seidel ou SGS (kfm12.eso (2D) et kfm14.eso (3D))
102 : C
103 : C************************************************************************
104 : C
105 : IMPLICIT INTEGER(I-N)
106 : -INC CCOPTIO
107 : -INC SMLMOTS
108 : -INC SMCHPOI
109 : POINTEUR MLMVIT.MLMOTS
110 : C
111 : C**** Variables de COOPTIO
112 : C
113 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
114 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
115 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
116 : C & ,IECHO, IIMPI, IOSPI
117 : C & ,IDIM
118 : C & ,MCOORD
119 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
120 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
121 : C & ,NORINC,NORVAL,NORIND,NORVAD
122 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI, IREFOR, ISAFOR
123 : CC
124 : INTEGER ICOEF, IDOMA, IRET, MELEF, MELEC, MELEFL, ICHPSU, INDIC
125 : & , NBCOMP, ICHPDI, IRN, IGN, IRETN, IGAMN, ICHPVO, ILIINC
126 : & , IDU, JGN, JGM, INORM, INEFMD, ICOND, MMODEL
127 : & , ICHLIM, MELLIM, ILIINP, NSOUPO, NJAC, IRES, IPROBL
128 : & , IVCO, MELTFA, ICOEV
129 : C
130 : REAL*8 DCFL, DTPS
131 : C
132 : CHARACTER*4 MOT
133 : CHARACTER*8 TYPE, LJACO(4)
134 : CHARACTER*(40) MESERR
135 : C
136 : DATA LJACO/'PJACO ','LJACOF ','LJACOB ','LJACOFB '/
137 : C
137 : CALL LIRMOT(LJACO,4,ICOEF,1)
138 : IF(IERR .NE. 0)GOTO 9999
139 : C
140 : C**********************************
141 : C**** Lecture de l'objet MODELE ***
142 : C**********************************
143 : C
144 : ICOND = 1
144 : CALL QUETYP(TYPE,ICOND,IRET)
145 :
146 : IF((IRET.EQ.0).AND.(TYPE.NE.'MMODEL'))THEN
147 : MOTERR(1:40)='MMODEL '
147 : CALL ERREUR(471)
148 : GOTO 9999
149 : ENDIF
149 : CALL LIROBJ('MMODEL',MMODEL,ICOND,IRET)
150 : IF(IERR.NE.0)GOTO 9999
150 : CALL LEKMOD(MMODEL,IDOMA,INEFMD)
151 : IF(IERR.NE.0)GOTO 9999
152 : C
153 : C**** CENTRE, FACE, FACEL, ELTFA
154 : C
154 : CALL LEKTAB(IDOMA,'CENTRE',MELEC)
155 : IF(IERR .NE. 0) GOTO 9999
156 : C
156 : CALL LEKTAB(IDOMA,'FACE',MELEF)
157 : IF(IERR .NE. 0) GOTO 9999
158 : C
158 : CALL LEKTAB(IDOMA,'FACEL',MELEFL)
159 : IF(IERR .NE. 0) GOTO 9999
160 : C
160 : CALL LEKTAB(IDOMA,'ELTFA',MELTFA)
161 : IF(IERR .NE. 0) GOTO 9999
162 : C
163 : C**** Lecture du CHPOINT contenant les surfaces des faces.
164 : C
164 : CALL LEKTAB(IDOMA,'XXSURFAC',ICHPSU)
165 : IF(IERR .NE. 0) GOTO 9999
166 : INDIC = 1
167 : NBCOMP = 1
168 : MOT = 'SCAL'
168 : CALL QUEPOI(ICHPSU, MELEF, INDIC, NBCOMP, MOT)
169 : IF(IERR .NE. 0) GOTO 9999
170 : C
171 : C**** Lecture du CHPOINT contenant les diametres minimums.
172 : C
172 : CALL LEKTAB(IDOMA,'XXDIEMIN',ICHPDI)
173 : IF(IERR .NE. 0) GOTO 9999
174 : INDIC = 1
175 : NBCOMP = 1
176 : MOT = 'SCAL'
176 : CALL QUEPOI(ICHPDI, MELEC, INDIC, NBCOMP, MOT)
177 : IF(IERR .NE. 0) GOTO 9999
178 :
179 : C
180 : C**** Lecture du CHPOINT contenant les normales aux faces
181 : C
182 : IF(IDIM .EQ. 2)THEN
183 : C Que les normales
183 : CALL LEKTAB(IDOMA,'XXNORMAF',INORM)
184 : IF(IERR .NE. 0) GOTO 9999
185 : JGN = 4
186 : JGM = 2
187 : SEGINI MLMVIT
188 : MLMVIT.MOTS(1) = 'UX '
189 : MLMVIT.MOTS(2) = 'UY '
189 : CALL QUEPO1(INORM, MELEF, MLMVIT)
190 : SEGSUP MLMVIT
191 : C
192 : ELSE
193 : C Les normales et les tangentes
194 : TYPE = ' '
194 : CALL ACMO(IDOMA,'MATROT',TYPE,INORM)
195 : IF (TYPE .NE. 'CHPOINT ') THEN
195 : CALL MATRAN(IDOMA,INORM)
196 : IF(IERR .NE. 0) GOTO 9999
197 : ENDIF
198 : JGN = 4
199 : JGM = 9
200 : SEGINI MLMVIT
201 : MLMVIT.MOTS(1) = 'UX '
202 : MLMVIT.MOTS(2) = 'UY '
203 : MLMVIT.MOTS(3) = 'UZ '
204 : MLMVIT.MOTS(4) = 'RX '
205 : MLMVIT.MOTS(5) = 'RY '
206 : MLMVIT.MOTS(6) = 'RZ '
207 : MLMVIT.MOTS(7) = 'MX '
208 : MLMVIT.MOTS(8) = 'MY '
209 : MLMVIT.MOTS(9) = 'MZ '
209 : CALL QUEPO1(INORM, MELEF, MLMVIT)
210 : SEGSUP MLMVIT
211 : ENDIF
212 : C
213 : C**** Lecture du CHPOINT contenant les volumes
214 : C
214 : CALL LEKTAB(IDOMA,'XXVOLUM',ICHPVO)
215 : IF(IERR .NE. 0) GOTO 9999
216 : INDIC = 1
217 : NBCOMP = 1
218 : MOT = 'SCAL'
218 : CALL QUEPOI(ICHPVO, MELEC, INDIC, NBCOMP, MOT)
219 : IF(IERR .NE. 0) GOTO 9999
220 : C
221 : C********************************
222 : C**** Fin table domaine *********
223 : C********************************
224 : C
225 : C**** La list des inconnues
226 : C
227 : TYPE='LISTMOTS'
227 : CALL LIROBJ(TYPE,ILIINC,1,IRET)
228 : IF(IERR .NE. 0) GOTO 9999
229 : C
230 : C**** On va lire les pointeurs des CHPOINTs
231 : C Lecture du CHPOINT residu
232 : C
233 : TYPE='CHPOINT'
233 : CALL LIROBJ(TYPE,IRES,1,IRET)
234 : IF(IERR.NE.0) GOTO 9999
235 : C
236 : C**** Controle du CHPOINT
237 : C
237 : CALL QUEPO1(IRES, MELEC, ILIINC)
238 : IF(IERR .NE. 0) GOTO 9999
239 : C
240 : C**** On va lire les pointeurs des CHPOINTs
241 : C Lecture du CHPOINT centre densité
242 : C
243 : TYPE='CHPOINT'
243 : CALL LIROBJ(TYPE,IRN,1,IRET)
244 : IF(IERR.NE.0) GOTO 9999
245 : C
246 : C**** Controle du CHPOINT
247 : C
248 : INDIC = 1
249 : NBCOMP = 1
250 : MOT = 'SCAL'
250 : CALL QUEPOI(IRN, MELEC, INDIC, NBCOMP, MOT)
251 : IF(IERR .NE. 0) GOTO 9999
252 : C
253 : C**** Lecture du CHPOINT QDM
254 : C
255 : TYPE='CHPOINT'
255 : CALL LIROBJ(TYPE,IGN,1,IRET)
256 : IF(IERR.NE.0) GOTO 9999
257 : C
258 : C**** Control du CHPOINT
259 : C
260 : INDIC = 1
261 : NBCOMP = IDIM
262 : JGN = 4
263 : JGM = IDIM
264 : SEGINI MLMOTS
265 : MLMOTS.MOTS(1) = 'UX '
266 : MLMOTS.MOTS(2) = 'UY '
267 : IF(IDIM .EQ. 3) MLMOTS.MOTS(3) = 'UZ '
267 : CALL QUEPO1(IGN, MELEC, MLMOTS)
268 : IF(IERR .NE. 0) GOTO 9999
269 : SEGSUP MLMOTS
270 : C
271 : C**** Lecture du CHPOINT centre energie totale
272 : C
273 : TYPE='CHPOINT'
273 : CALL LIROBJ(TYPE,IRETN,1,IRET)
274 : IF(IERR.NE.0) GOTO 9999
275 : C
276 : C**** Controle du CHPOINT
277 : C
278 : INDIC = 1
279 : NBCOMP = 1
280 : MOT = 'SCAL'
280 : CALL QUEPOI(IRETN, MELEC, INDIC, NBCOMP, MOT)
281 : IF(IERR .NE. 0) GOTO 9999
282 : C
283 : C**** Lecture du CHPOINT centre gamma
284 : C
285 : TYPE='CHPOINT'
285 : CALL LIROBJ(TYPE,IGAMN,1,IRET)
286 : IF(IERR.NE.0) GOTO 9999
287 : C
288 : C**** Controle du CHPOINT
289 : C
290 : INDIC = 1
291 : NBCOMP = 1
292 : MOT = 'SCAL'
292 : CALL QUEPOI(IGAMN, MELEC, INDIC, NBCOMP, MOT)
293 : IF(IERR .NE. 0) GOTO 9999
294 : C
295 : C**** Lecture du CHPOINT centre cutoff speed
296 : C
297 : TYPE='CHPOINT'
297 : CALL LIROBJ(TYPE,IVCO,1,IRET)
298 : IF(IERR.NE.0) GOTO 9999
299 : C
300 : C**** Controle du CHPOINT
301 : C
302 : INDIC = 1
303 : NBCOMP = 1
304 : MOT = 'SCAL'
304 : CALL QUEPOI(IVCO, MELEC, INDIC, NBCOMP, MOT)
305 : IF(IERR .NE. 0) GOTO 9999
306 : C
307 : C**** Lecture de DTPS
308 : C
308 : CALL LIRREE(DTPS,1,IRET)
309 : IF(IERR .NE. 0) GOTO 9999
310 : C
311 : C**** Lecture de le double de la CFL pour le temps dual
312 : C
312 : CALL LIRREE(DCFL,1,IRET)
313 : IF(IERR .NE. 0) GOTO 9999
314 : C
315 : C**** Lecture de nombre d'iterations
316 : C
316 : CALL LIRENT(NJAC,1,IRET)
317 : IF(IERR .NE. 0) GOTO 9999
318 : C
319 : C**** Lecture de conditions aux limites
320 : C
321 : IRET=0
321 : CALL LIRCHA(MOT,0,IRET)
322 : IF(IERR .NE. 0) GOTO 9999
323 : IF(IRET .NE. 0)THEN
324 : IF(MOT .EQ. 'CLIM')THEN
325 : C
326 : TYPE='LISTMOTS'
326 : CALL LIROBJ(TYPE,ILIINP,1,IRET)
327 : IF(IERR .NE. 0) GOTO 9999
328 : MLMOT1=ILIINP
329 : C
330 : TYPE='CHPOINT'
330 : CALL LIROBJ(TYPE,ICHLIM,1,IRET)
331 : IF(IERR.NE.0) GOTO 9999
332 : C
333 : MCHPOI = ICHLIM
334 : SEGACT MCHPOI
335 : NSOUPO = MCHPOI.IPCHP(/1)
336 : IF(NSOUPO .EQ. 0) THEN
337 : ICHLIM=0
338 : MELLIM=0
339 : ELSEIF(NSOUPO .GT. 1)THEN
340 : MOTERR(1:8) = 'CHAMPOIN'
341 : C
342 : C**************** Message d'erreur standard
343 : C 132 2
344 : C On veut un objet %m1:8 élémentaire
345 : C
345 : CALL ERREUR(132)
346 : GOTO 9999
347 : ELSE
348 : MSOUPO=MCHPOI.IPCHP(1)
349 : SEGACT MSOUPO
350 : MELLIM=MSOUPO.IGEOC
351 : SEGDES MSOUPO
352 : SEGDES MCHPOI
352 : CALL QUEPO1(ICHLIM, MELLIM, MLMOT1)
353 : IF(IERR.NE.0) GOTO 9999
354 : ENDIF
355 : C
356 : ELSE
356 : CALL REFUS
357 : ENDIF
358 : ELSE
359 : ICHLIM=0
360 : MELLIM=0
361 : ENDIF
362 : C
363 : C**** Le coeff. pour le calcul du rayon spectral visc.
364 : C
365 : C
366 : TYPE='CHPOINT'
366 : CALL LIROBJ(TYPE,ICOEV,1,IRET)
367 : IF(IERR.NE.0) GOTO 9999
368 : C
369 : C**** Controle du CHPOINT
370 : C
371 : INDIC = 1
372 : NBCOMP = 1
373 : MOT = 'SCAL'
373 : CALL QUEPOI(ICOEV, MELEC, INDIC, NBCOMP, MOT)
374 : IF(IERR .NE. 0) GOTO 9999
375 : C
376 : C**** Calcul de DUN
377 : C
378 : TYPE='CHPOINT'
378 : CALL KRCHP1(TYPE, MELEC, IDU, ILIINC)
379 : C
380 : C*** Choix de la procédure de relaxation (PJ ou SGS)
381 : C
382 : IF(ICOEF .EQ. 1)THEN
383 : IF(IDIM .EQ. 2)THEN
384 : C
385 : C*** Procédure de relaxation PJ
386 : C
386 : CALL KFM11(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
387 : & ICHPSU,ICHPDI,ICHPVO,INORM,
388 : & MELEC,MELEF,MELEFL,DTPS,DCFL,
389 : & MELLIM,ICHLIM,NJAC,ICOEV,
390 : & IDU,IPROBL)
391 : IF(IERR .NE. 0)GOTO 9999
392 : ELSE
392 : CALL KFM13(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
393 : & ICHPSU,ICHPDI,ICHPVO,INORM,
394 : & MELEC,MELEF,MELEFL,DTPS,DCFL,
395 : & MELLIM,ICHLIM,NJAC,ICOEV,
396 : & IDU,IPROBL)
397 : IF(IERR .NE. 0)GOTO 9999
398 : c 251 2
399 : c Tentative d'utilisation d'une option non implémentée
399 : c CALL ERREUR(251)
400 : c GOTO 9999
401 : ENDIF
402 : ELSE
403 : IF(IDIM .EQ. 2)THEN
404 : C
405 : C*** Procédure de relaxation SGS
406 : C
406 : CALL KFM12(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
407 : & ICHPSU,ICHPDI,ICHPVO,INORM,
408 : & MELEC,MELEF,MELEFL,MELTFA,DTPS,DCFL,
409 : & MELLIM,ICHLIM,NJAC,ICOEF,ICOEV,
410 : & IDU,IPROBL)
411 : IF(IERR .NE. 0)GOTO 9999
412 : ELSE
412 : CALL KFM14(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
413 : & ICHPSU,ICHPDI,ICHPVO,INORM,
414 : & MELEC,MELEF,MELEFL,MELTFA,DTPS,DCFL,
415 : & MELLIM,ICHLIM,NJAC,ICOEF,ICOEV,
416 : & IDU,IPROBL)
417 : C 251 2
418 : C Tentative d'utilisation d'une option non implémentée
418 : C CALL ERREUR(251)
419 : C GOTO 9999
420 : ENDIF
421 : ENDIF
422 : C
422 : CALL ECRENT(IPROBL)
423 : TYPE='CHPOINT '
423 : CALL ECROBJ(TYPE,IDU)
424 : C
425 : 9999 CONTINUE
426 : RETURN
427 : END
428 :
429 :
430 :
431 :
© Cast3M 2003 - All rights reserved.
Disclaimer