Commit 6c67d9a5 authored by Series Laurent's avatar Series Laurent

Commit initial

parents
FC = gfortran
LD = gfortran
FCFLAGS = -w
LDFLAGS =
LIBS =
OBJDIR = ./obj
EXE = run/mr_burgers
OBJS = obj/mod_common.o \
obj/mod_structure.o \
obj/mod_recherche.o \
obj/mod_fonction.o \
obj/mod_arbre.o \
obj/mod_initial.o \
obj/mod_print.o \
obj/mod_euler.o \
obj/cell_avg.o \
obj/mr_burgers1D.o
all: $(OBJDIR) $(EXE)
$(OBJDIR):
@echo "Creation du repertoire $@"
@mkdir -p $@
@echo
$(EXE): $(OBJS)
$(LD) $(LDFLAGS) $^ -o $@ $(LIBS)
obj/%.o : src/%.f90
$(FC) $(FCFLAGS) -J ./obj -c $< -o $@
obj/%.o : src/%.f
$(FC) $(FCFLAGS) -c $< -o $@
obj/mod_structure.o : src/mod_structure.f90 \
obj/mod_common.o
obj/mod_recherche.o : src/mod_recherche.f90 \
obj/mod_common.o \
obj/mod_structure.o
obj/mod_fonction.o : src/mod_fonction.f90 \
src/liberation.f90 \
src/viscosite.f90 \
obj/mod_common.o \
obj/mod_structure.o
obj/mod_arbre.o : src/mod_arbre.f90 \
src/encoding_moy.f90 \
src/encoding_det.f90 \
src/seuillage.f90 \
src/graduation.f90 \
src/graduation_local.f90 \
src/elagage.f90 \
src/chainage_support.f90 \
src/liberation.f90 \
src/liste_feuille.f90 \
src/feuille_fictive.f90 \
src/feuille_fictive_visc.f90 \
src/valeur_feuille_fictive.f90 \
src/decoding_moy.f90 \
obj/mod_common.o \
obj/mod_structure.o \
obj/mod_fonction.o \
obj/mod_recherche.o
obj/mod_initial.o : src/mod_initial.f90 \
src/lect_don.f90 \
src/chainage.f90 \
src/support_racine.f90 \
src/init_sol.f90 \
src/lect_reprise.f90 \
src/recherche.f90 \
src/init_flags.f90 \
src/multimesh.f90 \
obj/mod_common.o \
obj/mod_structure.o \
obj/mod_arbre.o
obj/mod_print.o : src/mod_print.f90 \
src/print_leaves.f90 \
src/print_reprise.f90 \
obj/mod_common.o \
obj/mod_structure.o \
obj/mod_recherche.o \
obj/mod_fonction.o
obj/mod_euler.o : src/mod_euler.f90 \
src/pas_de_temps.f90 \
src/undt_update.f90 \
src/u0_update.f90 \
src/flux_euler.f90 \
src/flux_osmp7.f90 \
src/flux_visc.f90 \
src/osmp7.f90 \
src/integration_euler.f90 \
src/integration1_visqueux.f90 \
src/integration2_visqueux.f90 \
src/mise_a_jour.f90 \
src/residu.f90 \
obj/mod_common.o \
obj/mod_structure.o \
obj/mod_recherche.o \
obj/mod_fonction.o
clean:
rm -f obj/*.o obj/*.mod run/mr_burgers
-1., 1. :left and right limits of the domain
1. :left value
0.8 :CFL number
1000000 :number of time steps ==> nmax
9 :max grid level ==> niv_max
1 :min grid level ==> niv_min
1 :width of reconstruction polyn ==> order = 2s+1
0.0001 :Harten's threshold ==> epsilon
0. :initial recorded time
.1 :time step for recorded time
0.4999999 :final time ==> stop simulation
subroutine cal_indice(niv_loc, ind_loc, ind)
use mod_common
IMPLICIT NONE
INTEGER :: l
INTEGER :: niv_loc
INTEGER, DIMENSION(ndim) :: ind, ind_loc
! =========== gestion des indices pour une recherche dans l'arbre =========
! +++++Application des conditions aux limites pour
! un indice (ind_loc) a un niveau (niv_loc) donnes
ind(:) = ind_loc(:)
do l = 1, ndim
! .....conditions aux limites ==> symetrie
if( caltype_deb(l) <= 1 .OR. caltype_fin(l) <= 1 ) then
do while(ind_loc(l) <= 0 .or. &
ind_loc(l) > 2**niv_loc*nbre_racine(l) )
if(ind_loc(l) <= 0 ) ind(l) = 1 - ind_loc(l)
if(ind_loc(l) > 2**niv_loc*nbre_racine(l) ) &
ind(l) = 2**(niv_loc+1)*nbre_racine(l) - ind_loc(l) + 1
ind_loc(l) = ind(l)
enddo
! .....conditions aux limites ==> periodicite
elseif( caltype_deb(l) == 2 .OR. caltype_fin(l) == 2 ) then
do while(ind_loc(l) <= 0 .or. &
ind_loc(l) > 2**niv_loc*nbre_racine(l) )
if(ind_loc(l) <= 0 ) &
ind(l) = 2**niv_loc*nbre_racine(l) + ind_loc(l)
if(ind_loc(l) > 2**niv_loc ) &
ind(l) = ind_loc(l) - 2**niv_loc*nbre_racine(l)
ind_loc(l) = ind(l)
enddo
endif
enddo
end subroutine cal_indice
c =================================================
subroutine cellave2(xlow,ylow,dx,dy,wl)
c =================================================
implicit none
double precision xlow, ylow, dx, dy, wl
integer i, iv
external fss, zeroin, fdisc
double precision zeroin, fss, fdisc
logical fl(5),alll,allr
double precision x(10),y(10),xx(5),yy(5)
double precision xc0, yc0, xc1, yc1, ss, area
common/fsscorn/ xc0,yc0,xc1,yc1
c
c # compute wl, fraction of cell that lies in left state.
c # For initial data with two states ql and qr separated by a
c # discontinuity. The curve along which the discontinuity lies is
c # specified by the function fdisc, which should return a value that
c # is negative on the side where ql lies and positive on the qr side.
c
c # xlow,ylow is the coordinate of the lower left corner of the cell.
c # dx, dy are grid spacing in x and y.
c
xx(1) = xlow
xx(2) = xlow
xx(3) = xlow+dx
xx(4) = xlow+dx
xx(5) = xx(1)
yy(1) = ylow
yy(2) = ylow+dy
yy(3) = ylow+dy
yy(4) = ylow
yy(5) = yy(1)
alll = .true.
allr = .true.
c
do 20 i=1,4
fl(i) = fdisc(xx(i),yy(i)) .lt. 0.d0
alll = alll .and. fl(i)
allr = allr .and. (.not. fl(i))
20 continue
fl(5) = fl(1)
c
if (alll) then
wl = 1.d0
return
endif
if (allr) then
wl = 0.d0
return
endif
c
iv = 0
do 40 i=1,4
if (fl(i)) then
iv = iv+1
x(iv) = xx(i)
y(iv) = yy(i)
endif
if (fl(i).neqv.fl(i+1)) then
iv = iv+1
xc0 = xx(i)
yc0 = yy(i)
xc1 = xx(i+1)
yc1 = yy(i+1)
ss = zeroin(0.d0, 1.d0, fss, 1d-8)
c write(27,*) 'xc,yc,ss:',xc0,yc0,xc1,yc1,ss
x(iv) = xx(i) + ss*(xx(i+1)-xx(i))
y(iv) = yy(i) + ss*(yy(i+1)-yy(i))
endif
40 continue
c
c # compute area:
c
if (iv.eq.0) then
wl = 0.d0
return
endif
c
x(iv+1) = x(1)
y(iv+1) = y(1)
area = 0.d0
do 50 i=1,iv
area = area + .5d0*(y(i)+y(i+1))*(x(i+1)-x(i))
c write(27,*) ' x,y:',x(i),y(i)
50 continue
c
wl = area / (dx*dy)
c write(27,*) 'area,wl:',area,wl
c
return
end
c
c
c
c
c
c =================================================
double precision function fss(s)
c =================================================
implicit none
double precision xc0, yc0, xc1, yc1, x, y, s, fdisc
common/fsscorn/ xc0,yc0,xc1,yc1
c
c # compute fdisc at distance s between corners (xc0,yc0) and (xc1,yc1)
c
x = xc0 + s*(xc1-xc0)
y = yc0 + s*(yc1-yc0)
fss = fdisc(x,y)
return
end
c =================================================
double precision function zeroin(ax,bx,f,tol)
c =================================================
implicit none
external f
double precision ax, bx, f, tol, eps, tol1, a, b, fa, fb
double precision c, fc, d, e, xm, s, p, q, r
c
c a zero of the function f(x) is computed in the interval ax,bx .
c (Standard routine from netlib)
c
c input..
c
c ax left endpoint of initial interval
c bx right endpoint of initial interval
c f function subprogram which evaluates f(x) for any x in
c the interval ax,bx
c tol desired length of the interval of uncertainty of the
c final result ( .ge. 0.0)
c
c
c output..
c
c zeroin abcissa approximating a zero of f in the interval ax,bx
c
c
c it is assumed that f(ax) and f(bx) have opposite signs
c without a check. zeroin returns a zero x in the given interval
c ax,bx to within a tolerance 4*macheps*dabs(x) + tol, where macheps
c is the relative machine precision.
c this function subprogram is a slightly modified translation of
c the algol 60 procedure zero given in richard brent, algorithms for
c minimization without derivatives, prentice - hall, inc. (1973).
c
c
c
c compute eps, the relative machine precision
c
eps = 1.d0
10 eps = eps/2.d0
tol1 = 1.0 + eps
if (tol1 .gt. 1.d0) go to 10
c
c initialization
c
a = ax
b = bx
fa = f(a)
fb = f(b)
c
c begin step
c
20 c = a
fc = fa
d = b - a
e = d
30 if (dabs(fc) .ge. dabs(fb)) go to 40
a = b
b = c
c = a
fa = fb
fb = fc
fc = fa
c
c convergence test
c
40 tol1 = 2.d0*eps*abs(b) + 0.5d0*tol
xm = .5d0*(c - b)
if (abs(xm) .le. tol1) go to 90
if (fb .eq. 0.d0) go to 90
c
c is bisection necessary
c
if (abs(e) .lt. tol1) go to 70
if (abs(fa) .le. abs(fb)) go to 70
c
c is quadratic interpolation possible
c
if (a .ne. c) go to 50
c
c linear interpolation
c
s = fb/fa
p = 2.0*xm*s
q = 1.0 - s
go to 60
c
c inverse quadratic interpolation
c
50 q = fa/fc
r = fb/fc
s = fb/fa
p = s*(2.d0*xm*q*(q - r) - (b - a)*(r - 1.d0))
q = (q - 1.d0)*(r - 1.0)*(s - 1.d0)
c
c adjust signs
c
60 if (p .gt. 0.d0) q = -q
p = abs(p)
c
c is interpolation acceptable
c
if ((2.d0*p) .ge. (3.d0*xm*q - abs(tol1*q))) go to 70
if (p .ge. abs(0.5d0*e*q)) go to 70
e = d
d = p/q
go to 80
c
c bisection
c
70 d = xm
e = d
c
c complete step
c
80 a = b
fa = fb
if (abs(d) .gt. tol1) b = b + d
c # DAC : NO 'sign' function in g95 (11/28/2005)
c if (abs(d) .le. tol1) b = b + sign(tol1,xm)
if (abs(d) .le. tol1) then
if (xm .ge. 0) then
b = b + tol1
else
b = b - xm
endif
endif
fb = f(b)
if ((fb*(fc/abs(fc))) .gt. 0.d0) go to 20
go to 30
c
c done
c
90 zeroin = b
return
end
c =================================================
function fdisc(x,y)
c =================================================
implicit double precision (a-h,o-z)
parameter (x0=0.5, y0=0., r0=0.25)
c
c # for computing cell averages for initial data that has a
c # discontinuity along some curve. fdisc should be negative to the
c # left of the curve and positive to the right
c
c # circle of radius r0:
fdisc = (x-x0)**2 + (y-y0)**2 - r0**2
c
return
end
recursive subroutine chainage(maille, Profondeur, Hauteur, &
Position, Pere_maille)
use mod_common
use mod_structure
IMPLICIT NONE
INTEGER :: i, k, l, m, n
INTEGER, intent(in) :: Profondeur, Hauteur
INTEGER, DIMENSION(ndim), intent(in) :: Position
INTEGER, DIMENSION(ndim) :: Pos
TYPE(structure_maille), POINTER :: maille
TYPE(structure_maille), POINTER :: Pere_maille
if (.not.associated(maille)) then
! .....On alloue la maille
allocate(maille)
! +++++Valorisation de la position
maille%maille_niv=Hauteur
maille%maille_i(:)=Position(:)
! +++++Valorisation de seuil suivant le niveau
! .....norme L1
maille%eps=epsilon * ( 2.**(ndim*(maille%maille_niv - niv_max)) )
! .....norme L2
!! maille%eps=epsilon * ( 2.**(ndim*(maille%maille_niv - niv_max)/2.) )
! +++++Positionnement des variables logiques
! .....Les points appartenant a l'arbre seront initialises plus tard
maille%arbre = .false.
! .....Les feuilles sont initialisees reelles (non fictives)
maille%feuille = .true.
! +++++Valorisation du chainage vers le pere
maille%pere => Pere_maille
! +++++Initialisation des liens
do l = 1, 2**ndim
nullify(maille%fils(l)%ptr)
enddo
do l = 1,ndim
nullify(maille%support(l)%precedent)
nullify(maille%support(l)%suivant)
enddo
endif
if( Profondeur > Hauteur) then
do l = 1, 2**ndim
k = 0
n = 0
do m = ndim, 1, -1
k = k+1
i = ((l-1) - n)/2**(ndim-k)
n = n + i*2**(ndim-k)
Pos(m) = 2*Position(m)+i-1
enddo
call chainage(maille%fils(l)%ptr, Profondeur, Hauteur+1, &
Pos, maille)
enddo
endif
end subroutine chainage
recursive subroutine chainage_support(Pere_maille)
use mod_common
use mod_structure
IMPLICIT NONE
INTEGER :: inc, k, l, m, n
INTEGER, DIMENSION(ndim) :: ijk
TYPE(structure_maille), POINTER :: Pere_maille
TYPE(structure_maille), POINTER :: maille
TYPE(structure_maille), POINTER :: courant
!
! ========== Parcours de l'arbre a partir de la racine =============
!
if( associated(Pere_maille%fils(1)%ptr) ) then
! +++++Allocation du support des fils
do l = 1, 2**ndim
! .....Pour chaque fils
maille => Pere_maille%fils(l)%ptr
! .....position de la maille
k = 0
n = 0
do m = ndim, 1, -1
k = k+1
ijk(m) = ((l-1) - n)/2**(ndim-k)
n = n + ijk(m)*2**(ndim-k)
enddo
! .....Allocation et Valorisation du support pour la maille courante
do m = 1, ndim
inc = 2**(m-1)
! *****Pour le fils gauche
if( ijk(m) == 0 ) then
! .....le suivant est son frere droit
maille%support(m)%suivant => Pere_maille%fils(l+inc)%ptr
! .....le precedent est son cousin droit
if( associated(Pere_maille%support(m)%precedent) ) then
courant => Pere_maille%support(m)%precedent
if( associated(courant%fils(l+inc)%ptr) ) then
maille%support(m)%precedent => courant%fils(l+inc)%ptr
else
nullify(maille%support(m)%precedent)
endif
else
nullify(maille%support(m)%precedent)
endif
! *****Pour le fils droit
else
! .....le precedent est son frere gauche
maille%support(m)%precedent => Pere_maille%fils(l-inc)%ptr
! .....le suivant est son cousin gauche
if(associated(Pere_maille%support(m)%suivant) ) then
courant => Pere_maille%support(m)%suivant
if( associated(courant%fils(l-inc)%ptr) ) then
maille%support(m)%suivant => courant%fils(l-inc)%ptr
else
nullify(maille%support(m)%suivant)
endif
else
nullify(maille%support(m)%suivant)
endif
endif
enddo
enddo
do l = 1, 2**ndim
call chainage_support(Pere_maille%fils(l)%ptr)
enddo
endif
end subroutine chainage_support
subroutine creation(maille, reelle)
use mod_common
use mod_structure
IMPLICIT