Skip to content

Commit 2ca2e66

Browse files
fluidnumerics-joefluidnumerics-joe
committed
Add initial draft of linear euler 3d with radiation example
* Need to add no-normal-flow boundary conditions, whicch will require implementation of tangent and binormal vectors. This will also necessitate an addition to the differential geometry notes on how we compute these additional boundary vectors.
1 parent 4d8e6aa commit 2ca2e66

File tree

4 files changed

+467
-0
lines changed

4 files changed

+467
-0
lines changed

‎examples/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,4 +66,5 @@ add_fortran_tests (
6666
"linear_euler2d_planewave_reflection.f90"
6767
"linear_euler2d_planewave_propagation.f90"
6868
"linear_shallow_water2d_nonormalflow.f90"
69+
"linear_euler3d_spherical_soundwave_radiation.f90"
6970
)

‎src/SELF_LinearEuler3D_t.f90

Lines changed: 269 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,269 @@
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 BUsLESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22+
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARIsLG IN ANY WAY OUT OF THE USE OF
23+
! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24+
!
25+
! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
26+
27+
module self_LinearEuler3D_t
28+
!! This module defines a class that can be used to solve the Linear Euler
29+
!! equations in 3-D. The Linear Euler Equations, here, are the Euler equations
30+
!! linearized about a motionless background state.
31+
!!
32+
!! The conserved variables are
33+
34+
!! \begin{equation}
35+
!! \vec{s} = \begin{pmatrix}
36+
!! \rho \\
37+
!! u \\
38+
!! v \\
39+
!! w \\
40+
!! p
41+
!! \end{pmatrix}
42+
!! \end{equation}
43+
!!
44+
!! The conservative flux is
45+
!!
46+
!! \begin{equation}
47+
!! \overleftrightarrow{f} = \begin{pmatrix}
48+
!! \rho_0 u \hat{x} + \rho_0 v \hat{y} + \rho_0 w \hat{z} \\
49+
!! \frac{p}{\rho_0} \hat{x} \\
50+
!! \frac{p}{\rho_0} \hat{y} \\
51+
!! \frac{p}{\rho_0} \hat{z} \\
52+
!! c^2 \rho_0 ( u \hat{x} + v \hat{y} + w \hat{z})
53+
!! \end{pmatrix}
54+
!! \end{equation}
55+
!!
56+
!! and the source terms are null.
57+
!!
58+
59+
use self_model
60+
use self_dgmodel3D
61+
use self_mesh
62+
63+
implicit none
64+
65+
type,extends(dgmodel3D) :: LinearEuler3D_t
66+
! Add any additional attributes here that are specific to your model
67+
real(prec) :: rho0 = 1.0_prec ! Reference density
68+
real(prec) :: c = 1.0_prec ! Sound speed
69+
real(prec) :: g = 0.0_prec ! gravitational acceleration (y-direction only)
70+
71+
contains
72+
procedure :: SourceMethod => sourcemethod_LinearEuler3D_t
73+
procedure :: SetNumberOfVariables => SetNumberOfVariables_LinearEuler3D_t
74+
procedure :: SetMetadata => SetMetadata_LinearEuler3D_t
75+
procedure :: entropy_func => entropy_func_LinearEuler3D_t
76+
!procedure :: hbc3D_NoNormalFlow => hbc3D_NoNormalFlow_LinearEuler3D_t
77+
procedure :: flux3D => flux3D_LinearEuler3D_t
78+
procedure :: riemannflux3D => riemannflux3D_LinearEuler3D_t
79+
!procedure :: source3D => source3D_LinearEuler3D_t
80+
procedure :: SphericalSoundWave => SphericalSoundWave_LinearEuler3D_t
81+
82+
endtype LinearEuler3D_t
83+
84+
contains
85+
86+
subroutine SetNumberOfVariables_LinearEuler3D_t(this)
87+
implicit none
88+
class(LinearEuler3D_t),intent(inout) :: this
89+
90+
this%nvar = 5
91+
92+
endsubroutine SetNumberOfVariables_LinearEuler3D_t
93+
94+
subroutine SetMetadata_LinearEuler3D_t(this)
95+
implicit none
96+
class(LinearEuler3D_t),intent(inout) :: this
97+
98+
call this%solution%SetName(1,"rho") ! Density
99+
call this%solution%SetUnits(1,"kg⋅m⁻³")
100+
101+
call this%solution%SetName(2,"u") ! x-velocity component
102+
call this%solution%SetUnits(2,"m⋅s⁻¹")
103+
104+
call this%solution%SetName(3,"v") ! y-velocity component
105+
call this%solution%SetUnits(3,"m⋅s⁻¹")
106+
107+
call this%solution%SetName(4,"w") ! z-velocity component
108+
call this%solution%SetUnits(4,"m⋅s⁻¹")
109+
110+
call this%solution%SetName(5,"P") ! Pressure
111+
call this%solution%SetUnits(5,"kg⋅m⁻¹⋅s⁻²")
112+
113+
endsubroutine SetMetadata_LinearEuler3D_t
114+
115+
pure function entropy_func_LinearEuler3D_t(this,s) result(e)
116+
!! The entropy function is the sum of kinetic and internal energy
117+
!! For the linear model, this is
118+
!!
119+
!! \begin{equation}
120+
!! e = \frac{1}{2} \left( \rho_0*( u^2 + v^2 ) + \frac{P^2}{\rho_0 c^2} \right)
121+
class(LinearEuler3D_t),intent(in) :: this
122+
real(prec),intent(in) :: s(1:this%nvar)
123+
real(prec) :: e
124+
125+
e = 0.5_prec*this%rho0*(s(2)*s(2)+s(3)*(3)+s(4)*s(4))+ &
126+
0.5_prec*(s(5)*s(5)/(this%rho0*this%c*this%c))
127+
128+
endfunction entropy_func_LinearEuler3D_t
129+
130+
! pure function hbc3D_NoNormalFlow_LinearEuler3D_t(this,s,nhat) result(exts)
131+
! class(LinearEuler3D_t),intent(in) :: this
132+
! real(prec),intent(in) :: s(1:this%nvar)
133+
! real(prec),intent(in) :: nhat(1:2)
134+
! real(prec) :: exts(1:this%nvar)
135+
! ! Local
136+
! integer :: ivar
137+
138+
! exts(1) = s(1) ! density
139+
! exts(2) = (nhat(2)**2-nhat(1)**2)*s(2)-2.0_prec*nhat(1)*nhat(2)*s(3) ! u
140+
! exts(3) = (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! v
141+
! exts(4) = (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! w
142+
! exts(5) = s(4) ! p
143+
144+
! endfunction hbc3D_NoNormalFlow_LinearEuler3D_t
145+
subroutine sourcemethod_LinearEuler3D_t(this)
146+
implicit none
147+
class(LinearEuler3D_t),intent(inout) :: this
148+
149+
return
150+
151+
endsubroutine sourcemethod_LinearEuler3D_t
152+
153+
pure function flux3D_LinearEuler3D_t(this,s,dsdx) result(flux)
154+
class(LinearEuler3D_t),intent(in) :: this
155+
real(prec),intent(in) :: s(1:this%nvar)
156+
real(prec),intent(in) :: dsdx(1:this%nvar,1:3)
157+
real(prec) :: flux(1:this%nvar,1:3)
158+
159+
flux(1,1) = this%rho0*s(2) ! density, x flux ; rho0*u
160+
flux(1,2) = this%rho0*s(3) ! density, y flux ; rho0*v
161+
flux(1,3) = this%rho0*s(4) ! density, y flux ; rho0*w
162+
163+
flux(2,1) = s(5)/this%rho0 ! x-velocity, x flux; p/rho0
164+
flux(2,2) = 0.0_prec ! x-velocity, y flux; 0
165+
flux(2,3) = 0.0_prec ! x-velocity, z flux; 0
166+
167+
flux(3,1) = 0.0_prec ! y-velocity, x flux; 0
168+
flux(3,2) = s(5)/this%rho0 ! y-velocity, y flux; p/rho0
169+
flux(3,3) = 0.0_prec ! y-velocity, z flux; 0
170+
171+
flux(4,1) = 0.0_prec ! z-velocity, x flux; 0
172+
flux(4,2) = 0.0_prec ! z-velocity, y flux; 0
173+
flux(4,3) = s(5)/this%rho0 ! z-velocity, z flux; p/rho0
174+
175+
flux(5,1) = this%c*this%c*this%rho0*s(2) ! pressure, x flux : rho0*c^2*u
176+
flux(5,2) = this%c*this%c*this%rho0*s(3) ! pressure, y flux : rho0*c^2*v
177+
flux(5,3) = this%c*this%c*this%rho0*s(4) ! pressure, y flux : rho0*c^2*w
178+
179+
endfunction flux3D_LinearEuler3D_t
180+
181+
pure function riemannflux3D_LinearEuler3D_t(this,sL,sR,dsdx,nhat) result(flux)
182+
!! Uses a local lax-friedrich's upwind flux
183+
!! The max eigenvalue is taken as the sound speed
184+
class(LinearEuler3D_t),intent(in) :: this
185+
real(prec),intent(in) :: sL(1:this%nvar)
186+
real(prec),intent(in) :: sR(1:this%nvar)
187+
real(prec),intent(in) :: dsdx(1:this%nvar,1:3)
188+
real(prec),intent(in) :: nhat(1:3)
189+
real(prec) :: flux(1:this%nvar)
190+
! Local
191+
real(prec) :: fL(1:this%nvar)
192+
real(prec) :: fR(1:this%nvar)
193+
real(prec) :: u,v,w,p,c,rho0
194+
195+
u = sL(2)
196+
v = sL(3)
197+
w = sL(4)
198+
p = sL(5)
199+
rho0 = this%rho0
200+
c = this%c
201+
fL(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density
202+
fL(2) = p*nhat(1)/rho0 ! u
203+
fL(3) = p*nhat(2)/rho0 ! v
204+
fL(4) = p*nhat(3)/rho0 ! w
205+
fL(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure
206+
207+
u = sR(2)
208+
v = sR(3)
209+
w = sR(4)
210+
p = sR(5)
211+
fR(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density
212+
fR(2) = p*nhat(1)/rho0 ! u
213+
fR(3) = p*nhat(2)/rho0 ! v'
214+
fR(4) = p*nhat(3)/rho0 ! w
215+
fR(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure
216+
217+
flux(1:5) = 0.5_prec*(fL(1:5)+fR(1:5))+c*(sL(1:5)-sR(1:5))
218+
219+
endfunction riemannflux3D_LinearEuler3D_t
220+
221+
subroutine SphericalSoundWave_LinearEuler3D_t(this,rhoprime,Lr,x0,y0,z0)
222+
!! This subroutine sets the initial condition for a weak blast wave
223+
!! problem. The initial condition is given by
224+
!!
225+
!! \begin{equation}
226+
!! \begin{aligned}
227+
!! \rho &= \rho_0 + \rho' \exp\left( -\ln(2) \frac{(x-x_0)^2 + (y-y_0)^2}{L_r^2} \right)
228+
!! u &= 0 \\
229+
!! v &= 0 \\
230+
!! E &= \frac{P_0}{\gamma - 1} + E \exp\left( -\ln(2) \frac{(x-x_0)^2 + (y-y_0)^2}{L_e^2} \right)
231+
!! \end{aligned}
232+
!! \end{equation}
233+
!!
234+
implicit none
235+
class(LinearEuler3D_t),intent(inout) :: this
236+
real(prec),intent(in) :: rhoprime,Lr,x0,y0,z0
237+
! Local
238+
integer :: i,j,k,iEl
239+
real(prec) :: x,y,z,rho,r,E
240+
241+
print*,__FILE__," : Configuring weak blast wave initial condition. "
242+
print*,__FILE__," : rhoprime = ",rhoprime
243+
print*,__FILE__," : Lr = ",Lr
244+
print*,__FILE__," : x0 = ",x0
245+
print*,__FILE__," : y0 = ",y0
246+
print*,__FILE__," : y0 = ",z0
247+
248+
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
249+
k=1:this%solution%N+1,iel=1:this%mesh%nElem)
250+
x = this%geometry%x%interior(i,j,k,iEl,1,1)-x0
251+
y = this%geometry%x%interior(i,j,k,iEl,1,2)-y0
252+
z = this%geometry%x%interior(i,j,k,iEl,1,3)-z0
253+
r = sqrt(x**2+y**2+z**2)
254+
255+
rho = (rhoprime)*exp(-log(2.0_prec)*r**2/Lr**2)
256+
257+
this%solution%interior(i,j,k,iEl,1) = rho
258+
this%solution%interior(i,j,k,iEl,2) = 0.0_prec
259+
this%solution%interior(i,j,k,iEl,3) = 0.0_prec
260+
this%solution%interior(i,j,k,iEl,4) = 0.0_prec
261+
this%solution%interior(i,j,k,iEl,5) = rho*this%c*this%c
262+
263+
enddo
264+
265+
call this%ReportMetrics()
266+
267+
endsubroutine SphericalSoundWave_LinearEuler3D_t
268+
269+
endmodule self_LinearEuler3D_t

‎src/cpu/SELF_LinearEuler3D.f90

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
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+
module self_LinearEuler3D
28+
29+
use self_LinearEuler3D_t
30+
31+
implicit none
32+
33+
type,extends(LinearEuler3D_t) :: LinearEuler3D
34+
endtype LinearEuler3D
35+
36+
endmodule self_LinearEuler3D

0 commit comments

Comments
 (0)