Skip to content

Commit b5fb449

Browse files
YonatanGMYonatanGM
authored
Add Fortran elastic-tube-1d solvers (#606)
1 parent 5373c4d commit b5fb449

File tree

12 files changed

+687
-0
lines changed

12 files changed

+687
-0
lines changed

‎changelog-entries/606.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
- Added Fortran solvers for the elastic-tube-1d tutorial (using the intrinsic Fortran bindings).

‎elastic-tube-1d/fluid-fortran/CMakeLists.txt‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
cmake_minimum_required(VERSION 3.10)
2+
project(ElasticTube LANGUAGES Fortran C)
3+
4+
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
5+
set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE)
6+
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
7+
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
8+
endif()
9+
message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}")
10+
11+
# Add bounds check and warnings in debug build
12+
if(CMAKE_BUILD_TYPE STREQUAL "Debug" AND CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
13+
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wall -Wextra -fbounds-check")
14+
endif()
15+
16+
find_package(precice 3 REQUIRED CONFIG)
17+
find_package(LAPACK REQUIRED)
18+
19+
add_executable(FluidSolver
20+
src/FluidComputeSolution.f90
21+
src/utilities.f90
22+
src/FluidSolver.f90
23+
)
24+
25+
target_link_libraries(FluidSolver PRIVATE precice::precice)
26+
target_link_libraries(FluidSolver PRIVATE ${LAPACK_LIBRARIES} ${LAPACK_LINKER_FLAGS})
27+

‎elastic-tube-1d/fluid-fortran/clean.sh‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
#!/usr/bin/env sh
2+
set -e -u
3+
4+
. ../../tools/cleaning-tools.sh
5+
6+
rm -rvf ./output/*.vtk
7+
clean_precice_logs .
8+
clean_case_logs .

‎elastic-tube-1d/fluid-fortran/run.sh‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
#!/usr/bin/env bash
2+
set -e -u
3+
4+
. ../../tools/log.sh
5+
exec > >(tee --append "$LOGFILE") 2>&1
6+
7+
if [ ! -d build ]; then
8+
mkdir build
9+
cmake -S . -B build
10+
cmake --build build
11+
fi
12+
13+
mkdir -p output
14+
15+
./build/FluidSolver ../precice-config.xml
16+
17+
close_log

‎elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90‎

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
module FluidComputeSolution
2+
implicit none
3+
integer, parameter :: dp = kind(1.0d0)
4+
5+
contains
6+
7+
subroutine fluid_compute_solution(velocity_old, pressure_old, &
8+
crossSectionLength_old, crossSectionLength, t, N, kappa, tau, &
9+
velocity, pressure, info)
10+
11+
real(dp), intent(in) :: velocity_old(:), pressure_old(:)
12+
real(dp), intent(in) :: crossSectionLength_old(:), crossSectionLength(:)
13+
real(dp), intent(in) :: t
14+
integer, intent(in) :: N
15+
real(dp), intent(in) :: kappa, tau
16+
real(dp), intent(inout) :: velocity(:), pressure(:)
17+
integer, intent(out) :: info
18+
19+
! Local variables
20+
integer :: i, k
21+
real(dp), parameter :: PI = 3.141592653589793_dp
22+
real(dp), parameter :: e = 10000.0_dp
23+
real(dp), parameter :: c_mk2 = e / 2.0_dp * sqrt(PI)
24+
real(dp), parameter :: u0 = 10.0_dp, ampl = 3.0_dp, frequency = 10.0_dp, &
25+
t_shift = 0.0_dp
26+
real(dp), parameter :: tolerance = 1.0e-15_dp
27+
integer, parameter :: max_iterations = 50
28+
29+
real(dp) :: alpha, L, dx, velocity_in, tmp2, norm_1, norm_2, norm
30+
31+
! LAPACK Variables
32+
integer :: nlhs, nrhs
33+
real(dp), allocatable :: Res(:)
34+
real(dp), allocatable :: LHS(:, :)
35+
integer, allocatable :: ipiv(:)
36+
37+
nlhs = 2*N + 2
38+
nrhs = 1
39+
40+
! Allocate arrays
41+
allocate (Res(2*N + 2))
42+
allocate (LHS(2*N + 2, 2*N + 2))
43+
allocate (ipiv(nlhs))
44+
45+
velocity = velocity_old
46+
pressure = pressure_old
47+
48+
! Stabilization intensity
49+
alpha = 0.0 !(N * kappa * tau) / (N * tau + 1);
50+
L = 10.0
51+
dx = L/kappa !1.0 / (N * kappa);
52+
53+
! Output status from dgesv (0 = success, < 0 = invalid argument, > 0 = singular matrix)
54+
info = 0
55+
56+
! Nonlinear solver loop
57+
do k = 1, max_iterations
58+
! Initialize residual vector
59+
Res = 0.0
60+
61+
! Compute residuals
62+
do i = 2, N ! Adjusted for 1-based indexing
63+
! Momentum
64+
Res(i) = (velocity_old(i)*crossSectionLength_old(i) - velocity(i)*crossSectionLength(i))*dx/tau
65+
Res(i) = Res(i) + 0.25*(-crossSectionLength(i + 1)*velocity(i)*velocity(i + 1) - &
66+
crossSectionLength(i)*velocity(i)*velocity(i + 1))
67+
Res(i) = Res(i) + 0.25*(-crossSectionLength(i + 1)*velocity(i)**2 - &
68+
crossSectionLength(i)*velocity(i)**2 + &
69+
crossSectionLength(i)*velocity(i - 1)*velocity(i) + &
70+
crossSectionLength(i - 1)*velocity(i - 1)*velocity(i))
71+
Res(i) = Res(i) + 0.25*(crossSectionLength(i - 1)*velocity(i - 1)**2 + &
72+
crossSectionLength(i)*velocity(i - 1)**2)
73+
Res(i) = Res(i) + 0.25*(crossSectionLength(i - 1)*pressure(i - 1) + &
74+
crossSectionLength(i)*pressure(i - 1) - &
75+
crossSectionLength(i - 1)*pressure(i) + &
76+
crossSectionLength(i + 1)*pressure(i) - &
77+
crossSectionLength(i)*pressure(i + 1) - &
78+
crossSectionLength(i + 1)*pressure(i + 1))
79+
80+
! Continuity
81+
Res(i + N + 1) = (crossSectionLength_old(i) - crossSectionLength(i))*dx/tau
82+
Res(i + N + 1) = Res(i + N + 1) + 0.25*(crossSectionLength(i - 1)*velocity(i - 1) + &
83+
crossSectionLength(i)*velocity(i - 1) + &
84+
crossSectionLength(i - 1)*velocity(i) - &
85+
crossSectionLength(i + 1)*velocity(i) - &
86+
crossSectionLength(i)*velocity(i + 1) - &
87+
crossSectionLength(i + 1)*velocity(i + 1))
88+
Res(i + N + 1) = Res(i + N + 1) + alpha*(pressure(i - 1) - 2.0*pressure(i) + pressure(i + 1))
89+
end do
90+
91+
! Boundary conditions
92+
velocity_in = u0 + ampl*sin(frequency*(t + t_shift)*PI)
93+
Res(1) = velocity_in - velocity(1)
94+
! Pressure Inlet is linearly interpolated
95+
Res(N + 2) = -pressure(1) + 2.0*pressure(2) - pressure(3)
96+
! Velocity Outlet is linearly interpolated
97+
Res(N + 1) = -velocity(N + 1) + 2.0*velocity(N) - velocity(N - 1)
98+
! Pressure Outlet is "non-reflecting"
99+
tmp2 = sqrt(c_mk2 - pressure_old(N + 1)/2.0) - &
100+
(velocity(N + 1) - velocity_old(N + 1))/4.0
101+
Res(2*N + 2) = -pressure(N + 1) + 2.0*(c_mk2 - tmp2**2)
102+
103+
! Compute residual norm
104+
norm_1 = sqrt(sum(Res**2))
105+
norm_2 = sqrt(sum(pressure**2) + sum(velocity**2))
106+
norm = norm_1/norm_2
107+
108+
if ((norm < tolerance .and. k > 1) .or. k > max_iterations) then
109+
exit
110+
end if
111+
112+
! Initialize the LHS matrix
113+
LHS = 0.0
114+
115+
! Populate LHS matrix
116+
do i = 2, N
117+
! Momentum, Velocity
118+
LHS(i, i - 1) = LHS(i, i - 1) + 0.25*(-2.0*crossSectionLength(i - 1)*velocity(i - 1) - &
119+
2.0*crossSectionLength(i)*velocity(i - 1) - &
120+
crossSectionLength(i)*velocity(i) - crossSectionLength(i - 1)*velocity(i))
121+
LHS(i, i) = LHS(i, i) + crossSectionLength(i)*dx/tau + &
122+
0.25*(crossSectionLength(i + 1)*velocity(i + 1) + &
123+
crossSectionLength(i)*velocity(i + 1) + &
124+
2.0*crossSectionLength(i + 1)*velocity(i) + &
125+
2.0*crossSectionLength(i)*velocity(i) - &
126+
crossSectionLength(i)*velocity(i - 1) - crossSectionLength(i - 1)*velocity(i - 1))
127+
LHS(i, i + 1) = LHS(i, i + 1) + 0.25*(crossSectionLength(i + 1)*velocity(i) + &
128+
crossSectionLength(i)*velocity(i))
129+
130+
! Momentum, Pressure
131+
LHS(i, N + 1 + i - 1) = LHS(i, N + 1 + i - 1) - 0.25*crossSectionLength(i - 1) - &
132+
0.25*crossSectionLength(i)
133+
LHS(i, N + 1 + i) = LHS(i, N + 1 + i) + 0.25*crossSectionLength(i - 1) - &
134+
0.25*crossSectionLength(i + 1)
135+
LHS(i, N + 1 + i + 1) = LHS(i, N + 1 + i + 1) + 0.25*crossSectionLength(i) + &
136+
0.25*crossSectionLength(i + 1)
137+
! Continuity, Velocity
138+
LHS(i + N + 1, i - 1) = LHS(i + N + 1, i - 1) - 0.25*crossSectionLength(i - 1) - &
139+
0.25*crossSectionLength(i)
140+
LHS(i + N + 1, i) = LHS(i + N + 1, i) - 0.25*crossSectionLength(i - 1) + &
141+
0.25*crossSectionLength(i + 1)
142+
LHS(i + N + 1, i + 1) = LHS(i + N + 1, i + 1) + 0.25*crossSectionLength(i) + &
143+
0.25*crossSectionLength(i + 1)
144+
145+
! Continuity, Pressure
146+
LHS(i + N + 1, N + 1 + i - 1) = LHS(i + N + 1, N + 1 + i - 1) - alpha
147+
LHS(i + N + 1, N + 1 + i) = LHS(i + N + 1, N + 1 + i) + 2.0*alpha
148+
LHS(i + N + 1, N + 1 + i + 1) = LHS(i + N + 1, N + 1 + i + 1) - alpha
149+
end do
150+
151+
! Boundary conditions in LHS
152+
! Velocity Inlet is prescribed
153+
LHS(1, 1) = 1.0
154+
! Pressure Inlet is linearly interpolated
155+
LHS(N + 2, N + 2) = 1.0
156+
LHS(N + 2, N + 3) = -2.0
157+
LHS(N + 2, N + 4) = 1.0
158+
! Velocity Outlet is linearly interpolated
159+
LHS(N + 1, N + 1) = 1.0
160+
LHS(N + 1, N) = -2.0
161+
LHS(N + 1, N - 1) = 1.0
162+
! Pressure Outlet is Non-Reflecting
163+
LHS(2*N + 2, 2*N + 2) = 1.0
164+
LHS(2*N + 2, N + 1) = -(sqrt(c_mk2 - pressure_old(N + 1)/2.0) - (velocity(N + 1) - velocity_old(N + 1))/4.0)
165+
166+
call dgesv(nlhs, nrhs, LHS, nlhs, ipiv, Res, nlhs, info)
167+
if (info /= 0) then
168+
write(*, *) "Linear Solver not converged!, Info: ", info
169+
end if
170+
171+
! Update velocity and pressure
172+
do i = 1, N + 1
173+
velocity(i) = velocity(i) + Res(i)
174+
pressure(i) = pressure(i) + Res(i + N + 1)
175+
end do
176+
end do
177+
178+
! Deallocate arrays
179+
deallocate(Res, LHS, ipiv)
180+
181+
end subroutine fluid_compute_solution
182+
end module FluidComputeSolution

0 commit comments

Comments
 (0)