-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathmain.f90
150 lines (111 loc) · 3.04 KB
/
main.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
program hydro1d
use datatypes_module
use grid_module
use params_module
use probparams_module, only: init_probparams
use init_module
use variables_module
use output_module
use runtime_diag_module
use bcs_module
use dt_module
use interface_states_godunov_module
use interface_states_plm_module
use interface_states_ppm_module
use interface_states_ppm_temp_module
use riemann_module
use update_module
implicit none
type(grid_t) :: grid
type(gridvar_t) :: U
type(gridedgevar_t) :: U_l, U_r, fluxes, godunov_state
type(gridvar_t) :: g
integer :: n, ng
real (kind=dp_t) :: t, dt, dt_new
! parse the inputs file
call init_params()
call init_probparams()
! set the number of ghostcells
if (godunov_type == 0) then
ng = 1
else
ng = 4
endif
! build the grid and storage for grid variables, interface states,
! and fluxes
call build(grid, nx, ng, xmin, xmax, xlboundary, xrboundary, is_spherical)
call build(U, grid, ncons)
call build(U_l, grid, ncons)
call build(U_r, grid, ncons)
call build(g, grid, 1)
call build(fluxes, grid, ncons)
call build(godunov_state, grid, nprim)
! set the initial conditions
n = 0
t = 0.0_dp_t
call init_data(U)
! construct the gravitational acceleration at cell-centers
if (do_gravity) then
call gravity(U, g)
endif
! set the boundary conditions
call fillBCs(U, g)
call output(U, t, n)
! main evolution loop
do while (t < tmax)
! compute the timestep
if (fixed_dt > 0) then
dt = fixed_dt
else
call compute_dt(U, n, dt_new)
if (n > 0) then
dt = min(dt_change*dt, dt_new)
else
dt = dt_new
endif
end if
if (t + dt > tmax) dt = tmax - t
! construct the interface states
if (godunov_type == 0) then
call make_interface_states_godunov(U, U_l, U_r, dt)
else if (godunov_type == 1) then
call make_interface_states_plm(U, g, U_l, U_r, dt)
else if (godunov_type == 2) then
if (ppm_temp) then
call make_interface_states_ppm_temp(U, g, U_l, U_r, dt)
else
call make_interface_states_ppm(U, g, U_l, U_r, dt)
endif
endif
! compute the fluxes
call solve_riemann(U_l, U_r, fluxes, godunov_state)
! do the conservative update
call update(U, g, fluxes, godunov_state, dt)
! construct the gravitational acceleration at cell-centers
if (do_gravity) then
call gravity(U, g)
endif
! set the boundary conditions
call fillBCs(U, g)
t = t + dt
n = n + 1
! output (if necessary)
print *, n, t, dt
if (plot_dt > 0.0_dp_t .and. &
mod(t - dt, plot_dt) > mod(t, plot_dt)) then
call output(U, t, n)
endif
! do any runtime diagnostics
call run_diag(U, t, n)
enddo
! final output
call output(U, t, n)
! clean-up
call destroy(g)
call destroy(U)
call destroy(U_l)
call destroy(U_r)
call destroy(grid)
call destroy(fluxes)
call destroy(godunov_state)
end program hydro1d