Skip to content

Commit 7cddfde

Browse files
eric2003eric2003
committed
add eno
1 parent 99e24a0 commit 7cddfde

File tree

522 files changed

+51652
-0
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

522 files changed

+51652
-0
lines changed

‎example/1d-linear-convection/eno2/fortran/01/CMakeLists.txt‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
cmake_minimum_required(VERSION 3.31)
2+
3+
project ( testprj )
4+
5+
set ( PRJ_COMPILE_FEATURES )
6+
set ( PRJ_COMPILE_DEFINITIONS )
7+
8+
enable_language(Fortran)
9+
10+
list ( APPEND PRJ_COMPILE_FEATURES cxx_std_20 )
11+
12+
if ( MSVC )
13+
set_property( DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY VS_STARTUP_PROJECT ${PROJECT_NAME} )
14+
endif()
15+
16+
add_executable( ${PROJECT_NAME}
17+
enoburgers.f90
18+
)
19+
20+
target_compile_features ( ${PROJECT_NAME}
21+
PRIVATE
22+
${PRJ_COMPILE_FEATURES}
23+
)
24+
25+
target_compile_definitions ( ${PROJECT_NAME}
26+
PRIVATE
27+
${PRJ_COMPILE_DEFINITIONS}
28+
)

‎example/1d-linear-convection/eno2/fortran/01/README.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
cmake ../ -T fortran=ifx

‎example/1d-linear-convection/eno2/fortran/01/enoburgers.f90‎

Lines changed: 267 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,267 @@
1+
module global
2+
implicit none
3+
integer, parameter :: nx = 40
4+
integer, parameter :: ighost = 10
5+
integer, parameter :: iorder = 2
6+
integer, parameter :: isize = iorder * ( iorder + 1 )
7+
real(8), parameter :: pi = 3.14159265358979323846
8+
integer il(0:nx),ir(0:nx)
9+
real(8) :: coef(-1:iorder-1,0:iorder-1)
10+
real(8) :: dd(0:ighost-1, -ighost:nx+ighost)
11+
real(8) :: up1_2m(0:nx), up1_2p(0:nx), flux(0:nx)
12+
real(8) :: res(1:nx)
13+
real(8) :: dx, dt
14+
end module global
15+
16+
module mesh_module
17+
use global, only: nx, ighost
18+
implicit none
19+
real(8) :: x(-ighost:nx+ighost)
20+
endmodule mesh_module
21+
22+
module field_module
23+
use global, only: nx, ighost
24+
implicit none
25+
real(8) :: pu(-ighost:nx+ighost), un(-ighost:nx+ighost)
26+
endmodule field_module
27+
28+
subroutine residual(u)
29+
use global
30+
implicit none
31+
real(8) :: u(-ighost:nx+ighost)
32+
integer i
33+
34+
call reconstruction(u)
35+
call engquist_osher_flux(up1_2m,up1_2p,flux,nx)
36+
do i = 1, nx
37+
res(i) = - ( flux(i) - flux(i-1) ) / dx
38+
enddo
39+
40+
end subroutine residual
41+
42+
subroutine reconstruction(u)
43+
use global
44+
implicit none
45+
real(8) :: u(-ighost:nx+ighost)
46+
integer :: i, j, m, k1, k2, l1, l2
47+
48+
!chose the stencil by ENO method
49+
do j = -ighost, nx + ighost
50+
dd(0,j)=u(j)
51+
enddo
52+
do i=1,iorder-1
53+
do j=-ighost,nx+ighost-1
54+
dd(i,j)=dd(i-1,j+1)-dd(i-1,j)
55+
enddo
56+
enddo
57+
58+
do j = 0, nx
59+
il(j) = j
60+
ir(j) = j + 1
61+
do i=1,iorder-1
62+
if( abs(dd(i,il(j)-1)) <= abs(dd(i,il(j))) ) then
63+
il(j)=il(j)-1
64+
endif
65+
if( abs(dd(i,ir(j)-1)) <= abs(dd(i,ir(j))) ) then
66+
ir(j)=ir(j)-1
67+
endif
68+
enddo
69+
enddo
70+
! reconstruction u(j+1_2)
71+
do j=0,nx
72+
k1=il(j)
73+
k2=ir(j)
74+
l1=j-k1
75+
l2=j-k2
76+
up1_2m(j)=0
77+
up1_2p(j)=0
78+
do m=0,iorder-1
79+
up1_2m(j)=up1_2m(j)+u(k1+m)*coef(l1,m)
80+
up1_2p(j)=up1_2p(j)+u(k2+m)*coef(l2,m)
81+
enddo
82+
enddo
83+
end subroutine reconstruction
84+
85+
!calculate numerical flux
86+
subroutine engquist_osher_flux(up1_2m,up1_2p,flux,nx)
87+
implicit none
88+
real(8) :: up1_2m(0:nx),up1_2p(0:nx),flux(0:nx)
89+
integer :: i, nx
90+
91+
do i = 0, nx
92+
if ( up1_2m(i) >= 0 ) then
93+
if ( up1_2p(i) >= 0 ) then
94+
flux(i) = 0.5 * up1_2m(i) * up1_2m(i)
95+
else
96+
flux(i) = 0.5 * ( up1_2m(i) * up1_2m(i) + up1_2p(i) * up1_2p(i) )
97+
endif
98+
else
99+
if ( up1_2p(i) >= 0 ) then
100+
flux(i) = 0
101+
else
102+
flux(i) = 0.5 * up1_2p(i) * up1_2p(i)
103+
endif
104+
endif
105+
enddo
106+
end subroutine engquist_osher_flux
107+
108+
subroutine boundary( u, nx, ighost )
109+
implicit none
110+
integer :: nx, ighost
111+
real(8) :: u(-ighost:nx+ighost)
112+
integer :: i
113+
114+
do i = 0, - ighost, - 1
115+
u( i ) = u( i + nx )
116+
enddo
117+
118+
do i = nx + 1, nx + ighost
119+
u( i ) = u( i - nx )
120+
enddo
121+
end subroutine boundary
122+
123+
subroutine update_oldfield(un, pu, nx, ighost)
124+
implicit none
125+
integer :: nx, ighost
126+
real(8) :: un(-ighost:nx+ighost), pu(-ighost:nx+ighost)
127+
integer :: i
128+
129+
do i = -ighost, nx + ighost
130+
un( i ) = pu( i )
131+
enddo
132+
end subroutine update_oldfield
133+
134+
subroutine init_coef
135+
use global
136+
implicit none
137+
real(8) :: values(isize) = [1.5d0, -0.5d0, 0.5d0, 0.5d0, -0.5d0, 1.5d0]
138+
integer :: i, j, icount
139+
140+
icount = 1
141+
do i = -1, iorder-1
142+
do j = 0, iorder-1
143+
coef(i, j) = values(icount)
144+
icount = icount + 1
145+
end do
146+
end do
147+
148+
end subroutine init_coef
149+
150+
subroutine init_mesh
151+
use global
152+
use mesh_module
153+
implicit none
154+
integer :: i
155+
156+
do i = -ighost, nx + ighost
157+
x(i) = ( i - 1 ) * dx + dx/2 - 1.0
158+
enddo
159+
160+
end subroutine init_mesh
161+
162+
subroutine init_field()
163+
use global
164+
use mesh_module
165+
use field_module
166+
implicit none
167+
integer :: i
168+
169+
do i = 1, nx
170+
pu(i) = 0.25 + 0.5 * sin( pi * x(i) )
171+
enddo
172+
173+
call boundary( pu, nx, ighost )
174+
call update_oldfield(un, pu, nx, ighost)
175+
176+
end subroutine init_field
177+
178+
subroutine runge_kutta_3()
179+
use global
180+
use field_module
181+
implicit none
182+
integer :: i
183+
real(8) :: c1, c2, c3
184+
185+
call residual(pu)
186+
do i = 1, nx
187+
pu(i) = pu(i) + dt * res(i)
188+
enddo
189+
call boundary( pu, nx, ighost )
190+
191+
call residual(pu)
192+
193+
do i = 1, nx
194+
pu(i) = 0.75 * un(i) + 0.25 * pu(i) + 0.25 * dt * res(i)
195+
enddo
196+
197+
call boundary( pu, nx, ighost )
198+
199+
call residual(pu)
200+
201+
c1 = 1.0 / 3.0
202+
c2 = 2.0 / 3.0
203+
c3 = 2.0 / 3.0
204+
205+
do i = 1, nx
206+
pu(i) = c1 * un(i) + c2 * pu(i) + c3 * dt * res(i)
207+
enddo
208+
call boundary( pu, nx, ighost )
209+
210+
call update_oldfield(un, pu, nx, ighost)
211+
212+
end subroutine runge_kutta_3
213+
214+
subroutine visualize
215+
use global
216+
use mesh_module
217+
use field_module
218+
implicit none
219+
integer :: i
220+
open(1,file='solution_total.plt',status='unknown')
221+
do i = -ighost, nx + ighost
222+
write(1,101) x(i),pu(i)
223+
enddo
224+
close(1)
225+
226+
open(2,file='solution.plt',status='unknown')
227+
do i = 1, nx
228+
write(2,101) x(i),pu(i)
229+
enddo
230+
close(2)
231+
101 format(1x,e20.10,e20.10)
232+
end subroutine visualize
233+
234+
program main
235+
use global
236+
use mesh_module
237+
use field_module
238+
implicit none
239+
integer :: i
240+
real(8) :: t, simu_time, xlen
241+
xlen = 2.0
242+
dx = xlen / nx
243+
dt = dx * 0.5
244+
write(*,*) 'nx=',nx
245+
write(*,*) 'dx=',dx
246+
write(*,*) 'Input T:'
247+
read(*,*) simu_time
248+
249+
call init_coef()
250+
call init_mesh()
251+
call init_field()
252+
253+
t = 0
254+
do while( t < simu_time )
255+
call runge_kutta_3()
256+
257+
t = t + dt
258+
if ( t + dt > simu_time ) then
259+
dt = simu_time - t
260+
endif
261+
enddo
262+
263+
write(*,*) t
264+
265+
call visualize()
266+
267+
end program main

0 commit comments

Comments
 (0)