* fichier : test.dgibi
******************************NONLINEAIRE DECENTRE******************************
***************************** CONVECTIF PUR *****************************
*
* On resout une equation hyperbolique non lineaire à l'aide
* d'une methode volumes finis decentre d'ordre 1 (Schema
* le plus simple possible)
* Elle s'ecrit d theta / dt + div (f(theta) u) = 0
*
*******************************************************************************
OPTI ECHO 0;
CRIT1 = 1.D-6 ;
GRAPH = 'N' ;
'SAUT' 'PAGE' ;
*
*---------------------------------------------------------------------
*
'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
'OPTI' 'ISOV' 'SURF' ;
*'OPTI' 'TRAC' 'PS' ;
*
*
*=========
* MAILLAGE
*=========
*
*
*- Création des points supports du contour du domaine, et des droites
*- passant par les centres et les faces pour le post-traitement.
*
L = 100.D0 ;
LS2 = L / 2.D0 ;
H = 1.D0 ;
HS2 = H / 2.D0 ;
X0 = 0.D0 ;
X1 = X0 + L ;
Y0 = 0.D0 ;
Y1 = Y0 + H ;
INUMX = 200 ;
INUMY = 1 ;
INUM1 = INUMX - 1 ;
Y01 = Y0 + Y1 * 0.5D0 ;
DX = X1 - X0 / INUMX ;
DX1 = DX / 2.D0 ;
XG = X0 + DX1 ;
XD = X1 - DX1 ;
*
A1 = X0 Y0 ;
A3 = X1 Y0 ;
D1 = X0 Y1 ;
D3 = X1 Y1 ;
B1 = X0 Y01 ;
B3 = X1 Y01 ;
C1 = XG Y01 ;
C3 = XD Y01 ;
P6 = LS2 Y01 ;
*
*- Création des DROITES frontieres
*
DRBAS = A3 'DROI' INUMX A1 ;
DRGAU = A1 'DROI' INUMY D1 ;
DRHAU = D1 'DROI' INUMX D3 ;
DRDRO = D3 'DROI' INUMY A3 ;
PELIM = DX1 / (5. * INUMX) ;
*
*- Creation maillage GEOMETRIQUE
*
PTOT1 = 'DALL' DRBAS DRGAU DRHAU DRDRO ;
PTOT2 = 'ORIE' PTOT1 ;
*
*- Creation maillage HYBRIDE y compris sous-objets (cond. limites)
*
DRMID = B1 'DROI' INUMX B3 ;
DRMIC = C1 'DROI' INUM1 C3 ;
EXT1 = 'MANU' 'POI1' B1 ;
*
QFTOT= 'CHAN' PTOT2 QUAF ;
QFGAU= 'CHAN' DRGAU QUAF ;
QFDRO= 'CHAN' DRDRO QUAF ;
QFHAU= 'CHAN' DRHAU QUAF ;
QFBAS= 'CHAN' DRBAS QUAF ;
ELIM PELIM (QFTOT ET QFGAU ET QFDRO ET DRMID ET DRMIC ET EXT1 ET
QFHAU ET QFBAS) ;
*
*================
* INITIALISATIONS
*================
*
* ----------------
* = MODELISATION =
* ----------------
MODHYB = 'MODL' QFTOT 'DARCY' 'ISOTROPE' ;
MODDRO = 'MODL' QFDRO 'DARCY' 'ISOTROPE' ;
MODGAU = 'MODL' QFGAU 'DARCY' 'ISOTROPE' ;
MODHAU = 'MODL' QFHAU 'DARCY' 'ISOTROPE' ;
MODBAS = 'MODL' QFBAS 'DARCY' 'ISOTROPE' ;
CHYB1 = 'DOMA' MODHYB 'SURFACE' ;
CHYB2 = 'DOMA' MODHYB 'NORMALE' ;
*
* ---------------------
* = Donnees physiques =
* ---------------------
*
T0 = 0.D0 ;
T1 = 1.D0 ;
VX1 = 1.D0 ; VY1 = 0.D0 ;
*VK = 1.D-2 ;
*
SPEED = 'MANU' 'CHPO' (DOMA MODHYB 'FACE') 2
'VX' VX1 'VY' VY1 'NATURE' 'DISCRET' ;
MOT1 = 'MOTS' 'VX' 'VY' ;
MOT2 = 'MOTS' 'UX' 'UY' ;
FLU1 = 'PSCA' SPEED CHYB2 MOT1 MOT2 ;
FLU2 = CHYB1 * FLU1 ;
FLU3 = 'NOMC' 'FLUX' FLU2 ;
*
* -----------------------
* = Donnees transitoire =
* -----------------------
* TETA : Parametre de le theta-méthode
* TMAX : Temps final
* TSUP : Temps pour conditions aux limites
* DELTAT : Pas de temps
*
TETA = 1.00D0 ;
TMIN = 0.D0 ;
TMAX = 40.00D0 ;
TSUP = 1.2D0 * TMAX ;
DELTAT = 0.2D0 ;
********************************************************************************
* CONDITIONS INITIALES
********************************************************************************
FNE= MANU CHPO (DOMA MODHYB 'CENTRE') 1 'SCAL' 0.D0 ;
********************************************************************************
* CONDITIONS AUX LIMITES;
********************************************************************************
CLIM1= MANU CHPO (DOMA MODGAU 'CENTRE') 1 'SCAL' 1.D0 ;
CLIM = CLIM1;
VHYB = DOMA MODHYB 'VOLUME';
MCHYB = 'DOMA' MODHYB 'ORIENTAT' ;
CONC2 = (-1.) * FLU3;
NOBP = TMAX/DELTAT;
NOBP = ENTI NOBP;
********************************************************************************
* BOUCLE EN TEMPS
*******************************************************************************
REPE BOUCT NOBP ;
NUM = (FNE**3)/2.;
DEN = (((FNE**3)/2.) + (((1-FNE)**3)/3.));
FCT = NUM*(DEN**-1);
AUX TERSC1 = 'DIVU' MODHYB CONC2 MCHYB FCT CLIM ;
DD= KCHA MODHYB TERSC1 'CHAM' ;
FNE=FNE+((TERSC1* DELTAT)*(VHYB**-1)) ;
FF= KCHA MODHYB FNE 'CHAM' ;
FIN BOUCT ;
********************************************************************************
SOLNUM = 'EVOL' 'VERT' 'CHPO' fne 'SCAL' DRMIC ;
ttot = tmax ;
k0 = prog 0.626 PAS 0.001 1.;
l2 = prog 375 * 1.;
k = inse k0 1 0.0;
********************************************************************************
* Solution analytique
********************************************************************************
beta = 18.;
a = (k0 '**' 2) '*' beta;
b = (l2 '-' k0) '**'2;
c = a '*' b;
l3 = prog 375 * 2.;
d = k0 '**' 3;
e = (k0 '**' 2) '*' 6.;
gam = -6.;
f = ( k0 '*' gam) '+' l3;
g = d '+' e '+' f;
h = g '**' 2;
dtanW = c '/' h;
dtanW = inse dtanW 1 1.3986075;
xt = dtanW*ttot;
xt = inse xt 1 100.;
k = inse k 1 0.0;
solanal = evol rouge manu xt k;
XX = EXTR solnum abs;
HH = EXTR solnum scal;
SOLPa = IPOL XX xt k ;
**************************
*erreur en norme Linfinie*
**************************
ER0 = HH - SOLPa ;
ERM = (MAXI (ABS ER0));
SI (ERM > 0.25);
'ERRE' 5 ;
SINON
'ERRE' 0 ;
FINSI;
MESSAGE ' ERREUR EN NORME LINFINIE ORDRE 1 = ' (MAXI (ABS ER0)) ;
********************
*erreur en norme L2*
********************
ER2 = ER0 * ER0 ;
ER1 = RESULT ER2 ;
ER = ER1/inumx;
ER = ER**0.5 ;
SI (ER > 0.0185);
'ERRE' 5 ;
SINON
'ERRE' 0 ;
FINSI;
MESSAGE ' ERREUR EN NORME L2 ABSOLUE ORDRE 1 = ' ER ;
SOLP = SOLPa * SOLPa;
SOLP = RESULT SOLP;
SOLP = SOLP**0.5;
ER1 = ER1**0.5;
ERR = ER1/SOLP;
SI (ERR > 0.0342);
'ERRE' 5 ;
SINON
'ERRE' 0 ;
FINSI;
*****************************
*erreur en norme L2 RELATIVE*
*****************************
MESSAGE ' ERREUR EN NORME L2 RELATIVE ORDRE 1 = ' ERR ;
TV =0.0 ;
j=0 ;
'REPETER' VARTOT (INUMX-1) ;
j=j+1 ;
i=j+1 ;
Vj = 'EXTR' hh j ;
Vj1 = 'EXTR' hh i ;
V=Vj1-Vj ;
TV = TV '+' ('ABS' V) ;
'FIN' VARTOT ;
'MESS' ' Variation totale ORDRE 1:' TV ;
TV =0.0 ;
j=0 ;
'REPETER' VARTOT (INUMX-1) ;
j=j+1 ;
i=j+1 ;
Vj = 'EXTR' solpa j ;
Vj1 = 'EXTR' solpa i ;
V=Vj1-Vj ;
TV = TV '+' ('ABS' V) ;
'FIN' VARTOT ;
*'MESS' ' Variation totale ORDRE 1:' TV ;
'SI' ('NEG' GRAPH 'N') ;
AV = SOLNUM 'ET' SOLANAL;
'DESS' AV 'MIMA' ;
FINSI;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer