Skip to content

Commit 1c28953

Browse files
zingalezingale
committed
add the MOL solver
1 parent 6acd793 commit 1c28953

19 files changed

Lines changed: 1553 additions & 0 deletions

‎compressible/MOL/Fortran/GNUmakefile‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
# will put all the build files (*.o, etc.) in a separate directory
2+
odir := _build
3+
4+
sources := $(wildcard *.f90)
5+
objects := $(addprefix $(odir)/, $(sources:.f90=.o))
6+
7+
ALL: euler
8+
9+
# we assume gfortran
10+
FFLAGS := -J $(odir) -I $(odir) -g -fbacktrace -fbounds-check -Wuninitialized -Wunused -ffpe-trap=invalid -finit-real=snan
11+
12+
# automagically find the interdependencies
13+
$(odir)/deps: $(sources)
14+
@if [ ! -d $(odir) ]; then mkdir -p $(odir); fi
15+
./util/dep.py --prefix $(odir)/ --search_path . $(sources) > $(odir)/deps
16+
17+
include $(odir)/deps
18+
19+
20+
$(odir)/%.o: %.f90
21+
@if [ ! -d $(odir) ]; then mkdir -p $(odir); fi
22+
gfortran $(FFLAGS) -c -o $@ $<
23+
24+
euler: $(objects)
25+
gfortran -o euler $(objects)
26+
27+
28+
# targets for cleaning up
29+
clean:
30+
rm -f $(odir)/*.o $(odir)/*.mod $(odir)/deps
31+
32+
realclean: clean
33+
@if [ -d $(odir) ]; then rmdir $(odir); echo "removing $(odir)"; fi
34+
rm -f euler

‎compressible/MOL/Fortran/advection.f90‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
module advection_module
2+
3+
use riemann_module
4+
5+
implicit none
6+
7+
contains
8+
9+
subroutine get_advective_source(U, k)
10+
11+
use grid_module
12+
use runtime_params_module
13+
14+
implicit none
15+
16+
real(dp_t), intent(in) :: U(a_lo:a_hi, NVAR)
17+
real(dp_t), intent(inout) :: k(a_lo:a_hi, NVAR)
18+
19+
real(dp_t) :: q(a_lo:a_hi, NVAR), dq(a_lo:a_hi, NVAR)
20+
real(dp_t) :: flux(a_lo:a_hi, NVAR)
21+
22+
integer :: i, n
23+
real(dp_t) :: dc, dl, dr, da, test
24+
real(dp_t) :: q_l(NVAR), q_r(NVAR)
25+
26+
! convert to primitive variables
27+
q(:, QRHO) = U(:, URHO)
28+
q(:, QU) = U(:, UMX)/q(:, QRHO)
29+
q(:, QP) = (U(:, UENER) - 0.5_dp_t * q(:, QRHO) * q(:, QU)**2)*(gamma - 1.0_dp_t)
30+
31+
! compute slopes
32+
do n = 1, NVAR
33+
do i = ilo-1, ihi+1
34+
dc = 0.5_dp_t * (q(i+1, n) - q(i-1, n))
35+
dl = q(i+1, n) - q(i, n)
36+
dr = q(i, n) - q(i-1, n)
37+
38+
test = dl * dr
39+
40+
if (test > 0.0_dp_t) then
41+
if (abs(dl) < abs(dr)) then
42+
da = dl
43+
else
44+
da = dr
45+
endif
46+
else
47+
da = 0.0
48+
endif
49+
50+
dq(i, n) = da
51+
52+
enddo
53+
enddo
54+
55+
! compute interface states -- this is a loop over interfaces
56+
do i = ilo, ihi+1
57+
58+
do n = 1, NVAR
59+
q_l(n) = q(i-1, n) + 0.5_dp_t * dq(i-1, n)
60+
q_r(n) = q(i, n) - 0.5_dp_t * dq(i, n)
61+
enddo
62+
63+
! solve Riemann problem
64+
call riemann(gamma, &
65+
q_l(QRHO), q_l(QU), q_l(QP), &
66+
q_r(QRHO), q_r(QU), q_r(QP), &
67+
flux(i, URHO), flux(i, UMX), flux(i, UENER))
68+
enddo
69+
70+
! compute the advective source
71+
do i = ilo, ihi
72+
do n = 1, NVAR
73+
k(i, n) = -(flux(i+1, n) - flux(i, n))/dx
74+
enddo
75+
enddo
76+
77+
end subroutine get_advective_source
78+
79+
end module advection_module

‎compressible/MOL/Fortran/datatypes.f90‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
module datatypes_module
2+
3+
implicit none
4+
5+
integer, parameter :: dp_t = selected_real_kind(15,307)
6+
7+
real (kind=dp_t) :: pi = 3.141592653589793238462643383279502884_dp_t
8+
9+
end module datatypes_module

‎compressible/MOL/Fortran/euler.f90‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
program euler
2+
3+
use advection_module
4+
use grid_module
5+
use init_module
6+
use io_module
7+
use riemann_module
8+
use runtime_params_module
9+
use timestep_module
10+
11+
implicit none
12+
13+
real(dp_t), allocatable :: U_old(:,:)
14+
real(dp_t), allocatable :: U_new(:,:)
15+
16+
! for the MOL righthand sides
17+
real(dp_t), allocatable :: k(:,:,:)
18+
19+
real(dp_t) :: t, dt
20+
21+
integer :: n, istep
22+
integer, parameter :: ng = 2
23+
24+
! read in runtime parameters
25+
call get_runtime_parameters()
26+
27+
! create the grid
28+
call init_grid(ng)
29+
30+
allocate(U_old(a_lo:a_hi, nvar))
31+
allocate(U_new(a_lo:a_hi, nvar))
32+
33+
U_old(:,:) = 0.0_dp_t
34+
U_new(:,:) = 0.0_dp_t
35+
36+
allocate(k(a_lo:a_hi, nvar, MOL_STAGES))
37+
38+
k(:,:,:) = 0.0_dp_t
39+
40+
! initialize the problem
41+
call prob_init(U_old)
42+
43+
! evolution loop
44+
t = 0.0
45+
istep = 0
46+
47+
call output(t, istep, U_old)
48+
49+
do while (t < tmax)
50+
51+
! compute the timestep
52+
call compute_dt(U_old, dt)
53+
54+
if (t + dt > tmax) then
55+
dt = tmax - t
56+
endif
57+
58+
! fill ghost cells
59+
call fill_ghost(U_old)
60+
61+
! get advective RHS
62+
call get_advective_source(U_old, k(:,:,1))
63+
64+
! predict state to next RK stage -- we'll temporarily store this in the
65+
! "new slot"
66+
do n = 1, nvar
67+
U_new(:, n) = U_old(:, n) + 0.5_dp_t * dt * k(:, n, 1)
68+
enddo
69+
70+
! fill ghost cells
71+
call fill_ghost(U_new)
72+
73+
! get advective RHS
74+
call get_advective_source(U_new, k(:,:,2))
75+
76+
! do final update
77+
do n = 1, nvar
78+
U_new(:, n) = U_old(:, n) + dt * k(:, n, 2)
79+
!U_new(:, n) = U_old(:, n) + dt * k(:, n, 1)
80+
enddo
81+
82+
t = t + dt
83+
istep = istep + 1
84+
85+
U_old(:, :) = U_new(:, :)
86+
87+
print *, istep, t, dt
88+
89+
enddo
90+
91+
call output(t, istep, U_new)
92+
93+
end program euler

‎compressible/MOL/Fortran/grid.f90‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
module grid_module
2+
3+
use datatypes_module
4+
5+
implicit none
6+
7+
! for indexing the conservative state arrays
8+
integer, parameter :: URHO = 1
9+
integer, parameter :: UMX = 2
10+
integer, parameter :: UENER = 3
11+
12+
! for indexing the primitive variable arrays
13+
integer, parameter :: QRHO = 1
14+
integer, parameter :: QU = 2
15+
integer, parameter :: QP = 3
16+
17+
integer, parameter :: NVAR = 3
18+
19+
integer, parameter :: MOL_STAGES = 2
20+
21+
real(dp_t), allocatable :: x(:)
22+
23+
real(dp_t) :: dx
24+
25+
integer :: ilo, ihi
26+
integer :: qx
27+
28+
integer :: a_lo, a_hi
29+
30+
contains
31+
32+
subroutine init_grid(ng)
33+
34+
use runtime_params_module
35+
36+
implicit none
37+
38+
integer, intent(in) :: ng
39+
40+
integer :: i
41+
42+
qx = nx + 2*ng
43+
44+
ilo = 0
45+
ihi = nx-1
46+
47+
a_lo = -ng
48+
a_hi = nx+ng-1
49+
50+
dx = (xmax - xmin)/nx
51+
52+
allocate(x(a_lo:a_hi))
53+
54+
do i = ilo-ng, ihi+ng
55+
x(i) = (dble(i) + 0.5_dp_t)*dx + xmin
56+
enddo
57+
58+
end subroutine init_grid
59+
60+
61+
subroutine fill_ghost(U)
62+
! just do zero-gradient
63+
64+
implicit none
65+
66+
integer :: n
67+
68+
real(dp_t), intent(inout) :: U(a_lo:a_hi, NVAR)
69+
70+
do n = 1, NVAR
71+
U(:ilo-1,n) = U(ilo,n)
72+
enddo
73+
74+
do n = 1, NVAR
75+
U(ihi+1:,n) = U(ihi,n)
76+
enddo
77+
78+
end subroutine fill_ghost
79+
80+
end module grid_module

‎compressible/MOL/Fortran/init.f90‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
module init_module
2+
3+
use datatypes_module
4+
5+
implicit none
6+
7+
contains
8+
9+
subroutine prob_init(U)
10+
11+
use grid_module
12+
use runtime_params_module
13+
14+
implicit none
15+
16+
real(dp_t) :: U(a_lo:a_hi, NVAR)
17+
18+
integer :: i
19+
20+
do i = ilo, ihi
21+
if (x(i) < 0.5_dp_t) then
22+
! left half
23+
U(i,URHO) = rho_l
24+
U(i,UMX) = rho_l * u_l
25+
U(i,UENER) = p_l/(gamma - 1.0_dp_t) + 0.5_dp_t * rho_l * u_l**2
26+
27+
else
28+
! left half
29+
U(i,URHO) = rho_r
30+
U(i,UMX) = rho_r * u_r
31+
U(i,UENER) = p_r/(gamma - 1.0_dp_t) + 0.5_dp_t * rho_r * u_r**2
32+
33+
endif
34+
35+
enddo
36+
37+
end subroutine prob_init
38+
39+
end module init_module

‎compressible/MOL/Fortran/inputs.double_rare‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
&params
2+
3+
rho_l = 1.0
4+
u_l = -2.0
5+
p_l = 0.4
6+
7+
rho_r = 1.0
8+
u_r = 2.0
9+
p_r = 0.4
10+
11+
tmax = 0.15
12+
13+
cfl = 0.8
14+
15+
nx = 128
16+
17+
/

‎compressible/MOL/Fortran/inputs.slow_shock‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
&params
2+
3+
rho_l = 5.6698
4+
u_l = -1.4701
5+
p_l = 100.0
6+
7+
rho_r = 1.0
8+
u_r = -10.5
9+
p_r = 1.0
10+
11+
tmax = 1.0
12+
13+
cfl = 0.8
14+
15+
nx = 128
16+
17+
/

‎compressible/MOL/Fortran/inputs.sod‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
&params
2+
3+
nx = 128
4+
5+
/

‎compressible/MOL/Fortran/inputs.strong_shock‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
&params
2+
3+
rho_l = 1.0
4+
u_l = 0.0
5+
p_l = 1000.0
6+
7+
rho_r = 1.0
8+
u_r = 0.0
9+
p_r = 0.01
10+
11+
tmax = 0.012
12+
13+
cfl = 0.8
14+
15+
nx = 128
16+
17+
/

0 commit comments

Comments
 (0)