Skip to content

Commit ff68de8

Browse files
committed
Use 2D StructuredMesh in shallow water
1 parent 6a7f3c7 commit ff68de8

File tree

3 files changed

+132
-5
lines changed

3 files changed

+132
-5
lines changed

examples/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,4 +80,5 @@ endfunction ()
8080

8181
add_fortran_tests (
8282
"burgers1d_shock.f90"
83+
"shallow_water_2d.f90"
8384
)

examples/shallow_water_2d.f90

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
2+
!
3+
! Maintainers : support@fluidnumerics.com
4+
! Official Repository : https://github.com/FluidNumerics/self/
5+
!
6+
! Copyright © 2024 Fluid Numerics LLC
7+
!
8+
! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
9+
!
10+
! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
11+
!
12+
! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in
13+
! the documentation and/or other materials provided with the distribution.
14+
!
15+
! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from
16+
! this software without specific prior written permission.
17+
!
18+
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
19+
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
20+
! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21+
! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22+
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
23+
! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24+
!
25+
! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
26+
27+
program shallow_water_2d_constant
28+
use self_data
29+
use self_shallow_water_2d
30+
31+
implicit none
32+
character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'euler'
33+
integer,parameter :: nvar = 3
34+
35+
integer,parameter :: controlDegree = 7
36+
integer,parameter :: targetDegree = 16
37+
real(prec),parameter :: H = 1.0_prec ! uniform resting depth
38+
real(prec),parameter :: g = 9.8_prec ! acceleration due to gravity
39+
real(prec),parameter :: dt = 1.0_prec*10.0_prec**(-5) ! time-step size
40+
real(prec),parameter :: endtime = 0.2_prec
41+
real(prec),parameter :: iointerval = 0.1_prec
42+
real(prec) :: e0,ef ! Initial and final entropy
43+
type(shallow_water_2d) :: modelobj
44+
type(Lagrange),target :: interp
45+
type(Mesh2D),target :: mesh
46+
integer :: bcids(1:4)
47+
type(SEMQuad),target :: geometry
48+
character(LEN=255) :: WORKSPACE
49+
50+
! Set boundary conditions
51+
bcids(1:4) = [SELF_BC_PRESCRIBED, & ! South
52+
SELF_BC_PRESCRIBED, & ! East
53+
SELF_BC_PRESCRIBED, & ! North
54+
SELF_BC_PRESCRIBED] ! West
55+
56+
! Create a uniform block mesh
57+
call mesh % StructuredMesh(10,10,2,2,0.05_prec,0.05_prec,bcids)
58+
59+
! Create an interpolant
60+
call interp%Init(N=controlDegree, &
61+
controlNodeType=GAUSS, &
62+
M=targetDegree, &
63+
targetNodeType=UNIFORM)
64+
65+
! Generate geometry (metric terms) from the mesh elements
66+
call geometry%Init(interp,mesh%nElem)
67+
call geometry%GenerateFromMesh(mesh)
68+
69+
! Initialize the model
70+
call modelobj%Init(nvar,mesh,geometry)
71+
modelobj%gradient_enabled = .true.
72+
73+
! Set the resting surface height and gravity
74+
modelobj%H = H
75+
modelobj%g = g
76+
77+
! Set the initial condition
78+
call modelobj%solution%SetEquation(1,'f = 1.0')
79+
call modelobj%solution%SetEquation(2,'f = 1.0')
80+
call modelobj%solution%SetEquation(3,'f = 1.0')
81+
call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
82+
83+
call modelobj%CalculateTendency()
84+
print*,"min, max (interior)", &
85+
minval(modelobj%solution%interior), &
86+
maxval(modelobj%solution%interior)
87+
88+
call modelobj%CalculateEntropy()
89+
e0 = modelobj%entropy
90+
91+
!Write the initial condition
92+
call modelobj%WriteModel()
93+
call modelobj%WriteTecplot()
94+
call modelobj%IncrementIOCounter()
95+
96+
! Set the model's time integration method
97+
call modelobj%SetTimeIntegrator(integrator)
98+
99+
! forward step the model to `endtime` using a time step
100+
! of `dt` and outputing model data every `iointerval`
101+
call modelobj%ForwardStep(endtime,dt,iointerval)
102+
103+
print*,"min, max (interior)", &
104+
minval(modelobj%solution%interior), &
105+
maxval(modelobj%solution%interior)
106+
107+
ef = modelobj%entropy
108+
109+
if(ef > e0) then
110+
print*,"Error: Final absmax greater than initial absmax! ",e0,ef
111+
stop 1
112+
endif
113+
! Clean up
114+
call modelobj%free()
115+
call mesh%free()
116+
call geometry%free()
117+
call interp%free()
118+
119+
endprogram shallow_water_2d_constant

test/shallow_water_2d_constant.f90

Lines changed: 12 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -43,12 +43,18 @@ program shallow_water_2d_constant
4343
type(shallow_water_2d) :: modelobj
4444
type(Lagrange),target :: interp
4545
type(Mesh2D),target :: mesh
46+
integer :: bcids(1:4)
4647
type(SEMQuad),target :: geometry
4748
character(LEN=255) :: WORKSPACE
4849

50+
! Set boundary conditions
51+
bcids(1:4) = [SELF_BC_PRESCRIBED, & ! South
52+
SELF_BC_PRESCRIBED, & ! East
53+
SELF_BC_PRESCRIBED, & ! North
54+
SELF_BC_PRESCRIBED] ! West
55+
4956
! Create a uniform block mesh
50-
call get_environment_variable("WORKSPACE",WORKSPACE)
51-
call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Block2D/Block2D_mesh.h5")
57+
call mesh % StructuredMesh(10,10,2,2,0.05_prec,0.05_prec,bcids)
5258

5359
! Create an interpolant
5460
call interp%Init(N=controlDegree, &
@@ -70,6 +76,8 @@ program shallow_water_2d_constant
7076

7177
! Set the initial condition
7278
call modelobj%solution%SetEquation(1,'f = 1.0')
79+
call modelobj%solution%SetEquation(2,'f = 1.0')
80+
call modelobj%solution%SetEquation(3,'f = 1.0')
7381
call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
7482

7583
call modelobj%CalculateTendency()
@@ -78,7 +86,6 @@ program shallow_water_2d_constant
7886
maxval(modelobj%solution%interior)
7987

8088
call modelobj%CalculateEntropy()
81-
call modelobj%ReportEntropy()
8289
e0 = modelobj%entropy
8390

8491
!Write the initial condition
@@ -100,8 +107,8 @@ program shallow_water_2d_constant
100107
ef = modelobj%entropy
101108

102109
if(ef > e0) then
103-
print*,"Error: Final absmax greater than initial absmax! ",e0,ef
104-
stop 1
110+
print*,"Error: Final absmax greater than initial absmax! ",e0,ef
111+
stop 1
105112
endif
106113
! Clean up
107114
call modelobj%free()

0 commit comments

Comments
 (0)