From 5da397d8b4845df84086c56a50a69994e09620a9 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Sun, 5 Jan 2025 09:27:22 +0100 Subject: [PATCH 01/24] Create initial folders and files --- elastic-tube-1d/fluid-fortran/CMakeLists.txt | 0 elastic-tube-1d/fluid-fortran/clean.sh | 0 elastic-tube-1d/fluid-fortran/run.sh | 0 elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 | 0 elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 | 0 elastic-tube-1d/fluid-fortran/src/utilities.f90 | 0 elastic-tube-1d/solid-fortran/CMakeLists.txt | 0 elastic-tube-1d/solid-fortran/clean.sh | 0 elastic-tube-1d/solid-fortran/run.sh | 0 elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 | 0 elastic-tube-1d/solid-fortran/src/SolidSolver.f90 | 0 11 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 elastic-tube-1d/fluid-fortran/CMakeLists.txt create mode 100644 elastic-tube-1d/fluid-fortran/clean.sh create mode 100644 elastic-tube-1d/fluid-fortran/run.sh create mode 100644 elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 create mode 100644 elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 create mode 100644 elastic-tube-1d/fluid-fortran/src/utilities.f90 create mode 100644 elastic-tube-1d/solid-fortran/CMakeLists.txt create mode 100644 elastic-tube-1d/solid-fortran/clean.sh create mode 100644 elastic-tube-1d/solid-fortran/run.sh create mode 100644 elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 create mode 100644 elastic-tube-1d/solid-fortran/src/SolidSolver.f90 diff --git a/elastic-tube-1d/fluid-fortran/CMakeLists.txt b/elastic-tube-1d/fluid-fortran/CMakeLists.txt new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/fluid-fortran/clean.sh b/elastic-tube-1d/fluid-fortran/clean.sh new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/fluid-fortran/run.sh b/elastic-tube-1d/fluid-fortran/run.sh new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/fluid-fortran/src/utilities.f90 b/elastic-tube-1d/fluid-fortran/src/utilities.f90 new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/solid-fortran/CMakeLists.txt b/elastic-tube-1d/solid-fortran/CMakeLists.txt new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/solid-fortran/clean.sh b/elastic-tube-1d/solid-fortran/clean.sh new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/solid-fortran/run.sh b/elastic-tube-1d/solid-fortran/run.sh new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 b/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 new file mode 100644 index 000000000..e69de29bb diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 new file mode 100644 index 000000000..e69de29bb From 1043c743b5fd26c48c05d8e97a0b914f6a791a92 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Sun, 5 Jan 2025 10:45:28 +0100 Subject: [PATCH 02/24] Add scripts, test with minimal solvers --- elastic-tube-1d/fluid-fortran/CMakeLists.txt | 31 ++++++++++++ elastic-tube-1d/fluid-fortran/clean.sh | 5 ++ elastic-tube-1d/fluid-fortran/run.sh | 12 +++++ .../fluid-fortran/src/FluidSolver.f90 | 48 +++++++++++++++++++ elastic-tube-1d/solid-fortran/CMakeLists.txt | 21 ++++++++ elastic-tube-1d/solid-fortran/clean.sh | 4 ++ elastic-tube-1d/solid-fortran/run.sh | 12 +++++ .../solid-fortran/src/SolidSolver.f90 | 42 ++++++++++++++++ 8 files changed, 175 insertions(+) mode change 100644 => 100755 elastic-tube-1d/fluid-fortran/clean.sh mode change 100644 => 100755 elastic-tube-1d/fluid-fortran/run.sh mode change 100644 => 100755 elastic-tube-1d/solid-fortran/clean.sh mode change 100644 => 100755 elastic-tube-1d/solid-fortran/run.sh diff --git a/elastic-tube-1d/fluid-fortran/CMakeLists.txt b/elastic-tube-1d/fluid-fortran/CMakeLists.txt index e69de29bb..a13db4621 100644 --- a/elastic-tube-1d/fluid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/fluid-fortran/CMakeLists.txt @@ -0,0 +1,31 @@ +cmake_minimum_required(VERSION 3.10) +project(ElasticTube_Fluid_Fortran LANGUAGES C Fortran DESCRIPTION "preCICE Fortran Fluid Solver") + +# If no build type is specified, default to Debug: +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE) + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS + "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +endif() +message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}") + +# Find preCICE library with Fortran bindings +find_package(precice 3 REQUIRED CONFIG) + +# Find LAPACK for linear solves +# Uncomment if LAPACK is needed +# find_package(LAPACK REQUIRED) + +# Add the executable +add_executable(FluidSolver +# Uncomment these if additional source files are required +# src/FluidComputeSolution.f90 +# src/utilities.f90 + src/FluidSolver.f90 +) + +# Link preCICE library +target_link_libraries(FluidSolver PRIVATE precice::precice) + +# Additional linking for LAPACK if needed +# target_link_libraries(FluidSolver PRIVATE LAPACK::LAPACK) diff --git a/elastic-tube-1d/fluid-fortran/clean.sh b/elastic-tube-1d/fluid-fortran/clean.sh old mode 100644 new mode 100755 index e69de29bb..d29b40321 --- a/elastic-tube-1d/fluid-fortran/clean.sh +++ b/elastic-tube-1d/fluid-fortran/clean.sh @@ -0,0 +1,5 @@ +#!/usr/bin/env bash +set -e -u + +rm -rf build || true +rm -rf output || true diff --git a/elastic-tube-1d/fluid-fortran/run.sh b/elastic-tube-1d/fluid-fortran/run.sh old mode 100644 new mode 100755 index e69de29bb..b0c7199f8 --- a/elastic-tube-1d/fluid-fortran/run.sh +++ b/elastic-tube-1d/fluid-fortran/run.sh @@ -0,0 +1,12 @@ +#!/usr/bin/env bash +set -e -u + +if [ ! -d build ]; then + mkdir build + cd build + cmake .. + cmake --build . + cd .. +fi + +./build/FluidSolver ../precice-config.xml diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index e69de29bb..7983ddf69 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -0,0 +1,48 @@ +program FluidSolver + implicit none + + character(len=512) :: configFileName + character(len=*), parameter :: solverName = 'Fluid' + character(len=*), parameter :: meshName = 'Fluid-Nodes-Mesh' + + integer :: rank, commsize, ongoing, bool + double precision :: dt + + ! We “declare” the preCICE native Fortran routines as external: + external :: precicef_create + external :: precicef_initialize + external :: precicef_is_coupling_ongoing + external :: precicef_get_max_time_step_size + external :: precicef_advance + external :: precicef_finalize + + write(*,*) "Starting minimal Fluid Solver..." + + if (command_argument_count() /= 1) then + write(*,*) "Usage: FluidSolver " + stop + endif + call get_command_argument(1, configFileName) + + rank = 0 + commsize = 1 + + ! create participant + call precicef_create(solverName, configFileName, rank, commsize) + + ! initialize + call precicef_initialize() + + ! minimal coupling loop + call precicef_is_coupling_ongoing(ongoing) + do while(ongoing /= 0) + call precicef_get_max_time_step_size(dt) + call precicef_advance(dt) + call precicef_is_coupling_ongoing(ongoing) + end do + + ! finalize + call precicef_finalize() + + write(*,*) "Exit FluidSolver" +end program FluidSolver diff --git a/elastic-tube-1d/solid-fortran/CMakeLists.txt b/elastic-tube-1d/solid-fortran/CMakeLists.txt index e69de29bb..ba56c7c3d 100644 --- a/elastic-tube-1d/solid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/solid-fortran/CMakeLists.txt @@ -0,0 +1,21 @@ +cmake_minimum_required(VERSION 3.10) +project(ElasticTube_Solid_Fortran LANGUAGES C Fortran DESCRIPTION "preCICE Fortran Solid Solver") + + +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE) + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS + "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +endif() +message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}") + + +find_package(precice REQUIRED CONFIG) + + +add_executable(SolidSolver + src/SolidSolver.f90 +) + + +target_link_libraries(SolidSolver PRIVATE precice::precice) diff --git a/elastic-tube-1d/solid-fortran/clean.sh b/elastic-tube-1d/solid-fortran/clean.sh old mode 100644 new mode 100755 index e69de29bb..bec0a3f61 --- a/elastic-tube-1d/solid-fortran/clean.sh +++ b/elastic-tube-1d/solid-fortran/clean.sh @@ -0,0 +1,4 @@ +#!/usr/bin/env bash +set -e -u + +rm -rf build || true diff --git a/elastic-tube-1d/solid-fortran/run.sh b/elastic-tube-1d/solid-fortran/run.sh old mode 100644 new mode 100755 index e69de29bb..e362cf904 --- a/elastic-tube-1d/solid-fortran/run.sh +++ b/elastic-tube-1d/solid-fortran/run.sh @@ -0,0 +1,12 @@ +#!/usr/bin/env bash +set -e -u + +if [ ! -d build ]; then + mkdir build + cd build + cmake .. + cmake --build . + cd .. +fi + +./build/SolidSolver ../precice-config.xml diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index e69de29bb..8f4d4f23f 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -0,0 +1,42 @@ +program SolidSolver + implicit none + + character(len=512) :: configFileName + character(len=*), parameter :: solverName = 'Solid' + + integer :: rank, commsize, ongoing + double precision :: dt + + external :: precicef_create + external :: precicef_initialize + external :: precicef_is_coupling_ongoing + external :: precicef_get_max_time_step_size + external :: precicef_advance + external :: precicef_finalize + + write(*,*) "Starting minimal Solid Solver..." + + if(command_argument_count() /= 1) then + write(*,*) "Usage: SolidSolver " + stop + end if + call get_command_argument(1, configFileName) + + rank = 0 + commsize = 1 + + call precicef_create(solverName, configFileName, rank, commsize) + + call precicef_initialize() + + call precicef_is_coupling_ongoing(ongoing) + do while(ongoing /= 0) + call precicef_get_max_time_step_size(dt) + call precicef_advance(dt) + call precicef_is_coupling_ongoing(ongoing) + end do + + call precicef_finalize() + + write(*,*) "Exit SolidSolver" +end program SolidSolver From 215449aca9054f7e44103a1a353471f9795741e0 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 16:47:28 +0100 Subject: [PATCH 03/24] Implementation wip --- elastic-tube-1d/fluid-fortran/CMakeLists.txt | 94 +++++-- elastic-tube-1d/fluid-fortran/clean.sh | 9 +- elastic-tube-1d/fluid-fortran/run.sh | 11 +- .../src/FluidComputeSolution.f90 | 203 +++++++++++++++ .../fluid-fortran/src/FluidSolver.f90 | 231 +++++++++++++++--- elastic-tube-1d/solid-fortran/CMakeLists.txt | 8 +- elastic-tube-1d/solid-fortran/clean.sh | 5 +- elastic-tube-1d/solid-fortran/run.sh | 5 + .../src/SolidComputeSolution.f90 | 28 +++ .../solid-fortran/src/SolidSolver.f90 | 135 +++++++--- 10 files changed, 638 insertions(+), 91 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/CMakeLists.txt b/elastic-tube-1d/fluid-fortran/CMakeLists.txt index a13db4621..e14fecfc4 100644 --- a/elastic-tube-1d/fluid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/fluid-fortran/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required(VERSION 3.10) -project(ElasticTube_Fluid_Fortran LANGUAGES C Fortran DESCRIPTION "preCICE Fortran Fluid Solver") +project(ElasticTube LANGUAGES Fortran C) -# If no build type is specified, default to Debug: +# Set default build type to Debug if not specified if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE) set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS @@ -9,23 +9,89 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) endif() message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}") -# Find preCICE library with Fortran bindings +# Find required packages find_package(precice 3 REQUIRED CONFIG) +find_package(LAPACK REQUIRED) -# Find LAPACK for linear solves -# Uncomment if LAPACK is needed -# find_package(LAPACK REQUIRED) - -# Add the executable +# Add executable add_executable(FluidSolver -# Uncomment these if additional source files are required -# src/FluidComputeSolution.f90 -# src/utilities.f90 + src/FluidComputeSolution.f90 + src/utilities.f90 src/FluidSolver.f90 ) -# Link preCICE library +# Link libraries target_link_libraries(FluidSolver PRIVATE precice::precice) +target_link_libraries(FluidSolver PRIVATE ${LAPACK_LIBRARIES} ${LAPACK_LINKER_FLAGS}) + +# ============================ +# Add Compiler Flags for Debugging +# ============================ + +# Define a function to add compiler flags based on compiler ID +function(add_fortran_debug_flags target) + if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") + # GNU Fortran (gfortran) specific flags + target_compile_options(${target} PRIVATE + $<$:-fbounds-check> + $<$:-fbacktrace> + $<$:-fcheck=all> + $<$:-g> + $<$:-Wall> + $<$:-Wextra> + ) + elseif(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") + # Intel Fortran (ifort) specific flags + target_compile_options(${target} PRIVATE + $<$:-check bounds> + $<$:-g> + $<$:-warn all> + ) + elseif(CMAKE_Fortran_COMPILER_ID MATCHES "PGI" OR + CMAKE_Fortran_COMPILER_ID MATCHES "NAG") + # PGI, NAG, and other Fortran compilers + target_compile_options(${target} PRIVATE + $<$:-C> + $<$:-g> + $<$:-Wall> + ) + else() + message(WARNING "Compiler ${CMAKE_Fortran_COMPILER_ID} not explicitly supported for debug flags.") + # Generic debug flags + target_compile_options(${target} PRIVATE + $<$:-g> + $<$:-Wall> + ) + endif() +endfunction() + +# Apply the debug flags to FluidSolver +add_fortran_debug_flags(FluidSolver) + +# ============================ +# Optionally, Add Linker Flags (if needed) +# ============================ + +# For example, to link LAPACK, you might need to specify additional flags +# However, typically, find_package(LAPACK) handles this +# If you have specific linker flags, you can add them as follows: +# target_link_libraries(FluidSolver PRIVATE ${LAPACK_LIBRARIES} ${LAPACK_LINKER_FLAGS}) + +# ============================ +# Enable Fortran Module Generation (Optional) +# ============================ + +# If you have modules that need to be visible to other targets, set the appropriate properties +# For example: +# set_target_properties(FluidSolver PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/modules) + +# ============================ +# Set CMake Policies (Optional) +# ============================ + +# To ensure compatibility with newer CMake versions +# cmake_policy(SET CMP0074 NEW) -# Additional linking for LAPACK if needed -# target_link_libraries(FluidSolver PRIVATE LAPACK::LAPACK) +# ============================ +# End of CMakeLists.txt +# ============================ diff --git a/elastic-tube-1d/fluid-fortran/clean.sh b/elastic-tube-1d/fluid-fortran/clean.sh index d29b40321..8d5d37784 100755 --- a/elastic-tube-1d/fluid-fortran/clean.sh +++ b/elastic-tube-1d/fluid-fortran/clean.sh @@ -1,5 +1,8 @@ -#!/usr/bin/env bash +#!/usr/bin/env sh set -e -u -rm -rf build || true -rm -rf output || true +. ../../tools/cleaning-tools.sh + +rm -rvf ./output/*.vtk +clean_precice_logs . +clean_case_logs . diff --git a/elastic-tube-1d/fluid-fortran/run.sh b/elastic-tube-1d/fluid-fortran/run.sh index b0c7199f8..2a0de7786 100755 --- a/elastic-tube-1d/fluid-fortran/run.sh +++ b/elastic-tube-1d/fluid-fortran/run.sh @@ -1,12 +1,15 @@ #!/usr/bin/env bash set -e -u +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + if [ ! -d build ]; then mkdir build - cd build - cmake .. - cmake --build . - cd .. + cmake -S . -B build + cmake --build build fi ./build/FluidSolver ../precice-config.xml + +close_log diff --git a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 index e69de29bb..c2fd73996 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 @@ -0,0 +1,203 @@ +module FluidComputeSolution + use, intrinsic :: iso_c_binding, only: c_double, c_int + implicit none + private + public :: fluidComputeSolutionSerial + +contains + + subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & + crossSectionLength_old, crossSectionLength, t, N, kappa, tau, & + velocity, pressure, info) + implicit none + + ! Function arguments + real(8), intent(in) :: velocity_old(:), pressure_old(:) + real(8), intent(in) :: crossSectionLength_old(:), crossSectionLength(:) + real(8), intent(in) :: t + integer, intent(in) :: N + real(8), intent(in) :: kappa, tau + real(8), intent(inout) :: velocity(:), pressure(:) + integer, intent(out) :: info + + ! Local variables + integer :: i, k, j + real(8), parameter :: PI = 3.141592653589793 + real(8), parameter :: E = 10000.0 + real(8), parameter :: c_mk2 = E / 2.0 * sqrt(PI) + real(8), parameter :: u0 = 10.0, ampl = 3.0, frequency = 10.0, & + t_shift = 0.0 + real(8), parameter :: tolerance = 1.0e-15 + integer, parameter :: max_iterations = 50 + + real(8) :: alpha, L, dx, velocity_in, tmp2, norm_1, norm_2, norm + + ! LAPACK Variables + integer :: nlhs, nrhs + real(8), allocatable :: Res(:) + real(8), allocatable :: LHS(:,:) + integer, allocatable :: ipiv(:) + + nlhs = 2 * N + 2 + nrhs = 1 + ! Allocate arrays + allocate(Res(2 * N + 2)) + allocate(LHS(2 * N + 2, 2 * N + 2)) + allocate(ipiv(nlhs)) + + velocity = velocity_old + pressure = pressure_old + + ! Stabilization intensity + alpha = 0.0 !(N * kappa * tau) / (N * tau + 1); + L = 10.0 + dx = L / kappa !1.0 / (N * kappa); + + ! Initialize result variable + info = 0 + print *, "CrossSectionLength:" + write(*, '(10F12.6)') crossSectionLength + print *, "Alpha:", alpha + + ! Nonlinear solver loop + do k = 1, max_iterations + ! Initialize residual vector + Res = 0.0 + + ! Compute residuals + do i = 2, N ! Adjusted for 1-based indexing + ! Momentum + Res(i) = (velocity_old(i) * crossSectionLength_old(i) - velocity(i) * crossSectionLength(i)) * dx / tau + Res(i) = Res(i) + 0.25 * (-crossSectionLength(i+1) * velocity(i) * velocity(i+1) - & + crossSectionLength(i) * velocity(i) * velocity(i+1)) + Res(i) = Res(i) + 0.25 * (-crossSectionLength(i+1) * velocity(i)**2 - & + crossSectionLength(i) * velocity(i)**2 + & + crossSectionLength(i) * velocity(i-1) * velocity(i) + & + crossSectionLength(i-1) * velocity(i-1) * velocity(i)) + Res(i) = Res(i) + 0.25 * (crossSectionLength(i-1) * velocity(i-1)**2 + & + crossSectionLength(i) * velocity(i-1)**2) + Res(i) = Res(i) + 0.25 * (crossSectionLength(i-1) * pressure(i-1) + & + crossSectionLength(i) * pressure(i-1) - & + crossSectionLength(i-1) * pressure(i) + & + crossSectionLength(i+1) * pressure(i) - & + crossSectionLength(i) * pressure(i+1) - & + crossSectionLength(i+1) * pressure(i+1)) + + ! Continuity + Res(i+N+1) = (crossSectionLength_old(i) - crossSectionLength(i)) * dx / tau + Res(i+N+1) = Res(i+N+1) + 0.25 * (crossSectionLength(i-1) * velocity(i-1) + & + crossSectionLength(i) * velocity(i-1) + & + crossSectionLength(i-1) * velocity(i) - & + crossSectionLength(i+1) * velocity(i) - & + crossSectionLength(i) * velocity(i+1) - & + crossSectionLength(i+1) * velocity(i+1)) + Res(i+N+1) = Res(i+N+1) + alpha * (pressure(i-1) - 2.0 * pressure(i) + pressure(i+1)) + end do + + ! Boundary conditions + velocity_in = u0 + ampl * sin(frequency * (t + t_shift) * PI) + Res(1) = velocity_in - velocity(1) + ! Pressure Inlet is linearly interpolated + Res(N+2) = -pressure(1) + 2.0 * pressure(2) - pressure(3) + ! Velocity Outlet is linearly interpolated + Res(N+1) = -velocity(N+1) + 2.0 * velocity(N) - velocity(N-1) + ! Pressure Outlet is "non-reflecting" + tmp2 = sqrt(c_mk2 - pressure_old(N+1) / 2.0) - & + (velocity(N+1) - velocity_old(N+1)) / 4.0 + Res(2*N+2) = -pressure(N+1) + 2.0 * (c_mk2 - tmp2**2) + + ! Compute residual norm + norm_1 = sqrt(sum(Res**2)) + norm_2 = sqrt(sum(pressure**2) + sum(velocity**2)) + norm = norm_1 / norm_2 + + if ((norm < tolerance .and. k > 1) .or. k > max_iterations) then + exit + end if + + ! Initialize the LHS matrix + LHS = 0.0 + + ! Populate LHS matrix + do i = 2, N + ! Momentum, Velocity + LHS(i, i-1) = LHS(i, i-1) + 0.25 * (-2.0 * crossSectionLength(i-1) * velocity(i-1) - & + 2.0 * crossSectionLength(i) * velocity(i-1) - & + crossSectionLength(i) * velocity(i) - crossSectionLength(i-1) * velocity(i)) + LHS(i, i) = LHS(i, i) + crossSectionLength(i) * dx / tau + & + 0.25 * (crossSectionLength(i+1) * velocity(i+1) + & + crossSectionLength(i) * velocity(i+1) + & + 2.0 * crossSectionLength(i+1) * velocity(i) + & + 2.0 * crossSectionLength(i) * velocity(i) - & + crossSectionLength(i) * velocity(i-1) - crossSectionLength(i-1) * velocity(i-1)) + LHS(i, i+1) = LHS(i, i+1) + 0.25 * (crossSectionLength(i+1) * velocity(i) + & + crossSectionLength(i) * velocity(i)) + + ! Momentum, Pressure + LHS(i, N+1+i-1) = LHS(i, N+1+i-1) - 0.25 * crossSectionLength(i-1) - & + 0.25 * crossSectionLength(i) + LHS(i, N+1+i) = LHS(i, N+1+i) + 0.25 * crossSectionLength(i-1) - & + 0.25 * crossSectionLength(i+1) + LHS(i, N+1+i+1) = LHS(i, N+1+i+1) + 0.25 * crossSectionLength(i) + & + 0.25 * crossSectionLength(i+1) + ! Continuity, Velocity + LHS(i+N+1, i-1) = LHS(i+N+1, i-1) - 0.25 * crossSectionLength(i-1) - & + 0.25 * crossSectionLength(i) + LHS(i+N+1, i) = LHS(i+N+1, i) - 0.25 * crossSectionLength(i-1) + & + 0.25 * crossSectionLength(i+1) + LHS(i+N+1, i+1) = LHS(i+N+1, i+1) + 0.25 * crossSectionLength(i) + & + 0.25 * crossSectionLength(i+1) + + ! Continuity, Pressure + LHS(i+N+1, N+1+i-1) = LHS(i+N+1, N+1+i-1) - alpha + LHS(i+N+1, N+1+i) = LHS(i+N+1, N+1+i) + 2.0 * alpha + LHS(i+N+1, N+1+i+1) = LHS(i+N+1, N+1+i+1) - alpha + end do + + + ! Boundary conditions in LHS + ! Velocity Inlet is prescribed + LHS(1, 1) = 1.0 + ! Pressure Inlet is linearly interpolated + LHS(N+2, N+2) = 1.0 + LHS(N+2, N+3) = -2.0 + LHS(N+2, N+4) = 1.0 + ! Velocity Outlet is linearly interpolated + LHS(N+1, N+1) = 1.0 + LHS(N+1, N) = -2.0 + LHS(N+1, N-1) = 1.0 + ! Pressure Outlet is Non-Reflecting + LHS(2*N+2, 2*N+2) = 1.0 + LHS(2*N+2, N+1) = -(sqrt(c_mk2 - pressure_old(N+1) / 2.0) - (velocity(N+1) - velocity_old(N+1)) / 4.0) + + + ! Solve Ax = b using LAPACK + ! print *, "LHS Matrix (size: ", size(LHS, 1), "x", size(LHS, 2), "):" + ! do i = 1, size(LHS, 1) + ! write(*, '(10F12.6)') (LHS(i, j), j = 1, size(LHS, 2)) + ! end do + ! print *, "Diagonal elements of LHS:" + ! do i = 1, size(LHS, 1) + ! print *, "LHS(", i, ",", i, ") =", LHS(i, i) + ! end do + print *, "LHS Matrix (size: ", size(LHS, 1), "x", size(LHS, 2), "):" + do i = 1, size(LHS, 1) + write(*, '(10F8.2)') (LHS(i, j), j = 1, size(LHS, 2)) + end do + + call dgesv(nlhs, nrhs, LHS, nlhs, ipiv, Res, nlhs, info) + if (info /= 0) then + return + end if + + ! Update velocity and pressure + do i = 1, N+1 + velocity(i) = velocity(i) + Res(i) + pressure(i) = pressure(i) + Res(i+N+1) + end do + end do + + ! Deallocate arrays + deallocate(Res, LHS, ipiv) + end subroutine fluidComputeSolutionSerial +end module FluidComputeSolution \ No newline at end of file diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index 7983ddf69..e8c8c2768 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -1,48 +1,211 @@ -program FluidSolver - implicit none +PROGRAM FluidSolver + USE FluidComputeSolution, ONLY: fluidComputeSolutionSerial + USE Utilities_mod, ONLY: write_vtk + IMPLICIT NONE - character(len=512) :: configFileName - character(len=*), parameter :: solverName = 'Fluid' - character(len=*), parameter :: meshName = 'Fluid-Nodes-Mesh' + ! Variable Declarations + CHARACTER(LEN=512) :: configFileName + CHARACTER(LEN=50) :: solverName + CHARACTER(LEN=50) :: meshName, pressureName, crossSectionLengthName + CHARACTER(LEN=256) :: outputFilePrefix + INTEGER :: rank, commsize, ongoing, dimensions, bool + INTEGER :: domainSize, chunkLength + INTEGER :: i, j, info + DOUBLE PRECISION :: dt, t, cellwidth + DOUBLE PRECISION, ALLOCATABLE :: pressure(:), pressure_old(:) + DOUBLE PRECISION, ALLOCATABLE :: crossSectionLength(:), crossSectionLength_old(:) + DOUBLE PRECISION, ALLOCATABLE :: velocity(:), velocity_old(:) + INTEGER, ALLOCATABLE :: vertexIDs(:) + INTEGER :: out_counter + DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793d0 + DOUBLE PRECISION :: kappa, L + DOUBLE PRECISION :: r0, a0, u0, ampl, frequency, t_shift, p0, vel_in_0 + DOUBLE PRECISION, ALLOCATABLE :: grid(:) - integer :: rank, commsize, ongoing, bool - double precision :: dt + ! Start of Program + WRITE (*,*) 'Fluid: Starting Fortran solver...' - ! We “declare” the preCICE native Fortran routines as external: - external :: precicef_create - external :: precicef_initialize - external :: precicef_is_coupling_ongoing - external :: precicef_get_max_time_step_size - external :: precicef_advance - external :: precicef_finalize + ! Command-Line Argument Parsing + IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN + WRITE (*,*) "" + WRITE (*,*) "Fluid: Usage: FluidSolver " + STOP -1 + END IF - write(*,*) "Starting minimal Fluid Solver..." + CALL getarg(1, configFileName) - if (command_argument_count() /= 1) then - write(*,*) "Usage: FluidSolver " - stop - endif - call get_command_argument(1, configFileName) + solverName = 'Fluid' + outputFilePrefix = './output/out_fluid' + ! Initialize preCICE Interface rank = 0 commsize = 1 + CALL precicef_create(solverName, configFileName, rank, commsize) + WRITE (*,*) "preCICE configured..." - ! create participant - call precicef_create(solverName, configFileName, rank, commsize) + ! Define Mesh and Data Names + meshName = "Fluid-Nodes-Mesh" + pressureName = "Pressure" + crossSectionLengthName = "CrossSectionLength" - ! initialize - call precicef_initialize() + domainSize = 4 + chunkLength = domainSize + 1 + kappa = 100.0d0 + L = 10.0d0 - ! minimal coupling loop - call precicef_is_coupling_ongoing(ongoing) - do while(ongoing /= 0) - call precicef_get_max_time_step_size(dt) - call precicef_advance(dt) - call precicef_is_coupling_ongoing(ongoing) + ! Get Mesh Dimensions from preCICE + CALL precicef_get_mesh_dimensions(meshName, dimensions) + + ! Allocate Arrays + ALLOCATE(vertexIDs(chunkLength)) + ALLOCATE(pressure(chunkLength)) + ALLOCATE(pressure_old(chunkLength)) + ALLOCATE(crossSectionLength(chunkLength)) + ALLOCATE(crossSectionLength_old(chunkLength)) + ALLOCATE(velocity(chunkLength)) + ALLOCATE(velocity_old(chunkLength)) + ALLOCATE(grid(dimensions * chunkLength)) + + ! Initialize vertexIDs (0-based IDs) + DO i = 1, chunkLength + vertexIDs(i) = i - 1 + END DO + + ! Initialize Physical Parameters + r0 = 1.0d0 / SQRT(PI) + a0 = r0**2 * PI + u0 = 10.0d0 + ampl = 3.0d0 + frequency = 10.0d0 + t_shift = 0.0d0 + p0 = 0.0d0 + vel_in_0 = u0 + ampl * SIN(frequency * (t_shift) * PI) + + ! Initialize Data Arrays + pressure = p0 + pressure_old = pressure + crossSectionLength = a0 + crossSectionLength_old = crossSectionLength + velocity = vel_in_0 + velocity_old = velocity + + ! Initialize Grid Coordinates + cellwidth = L / REAL(domainSize, KIND=8) + DO i = 1, chunkLength + DO j = 1, dimensions + IF (j == 1) THEN + grid((i - 1) * dimensions + j) = REAL(i - 1, KIND=8) * cellwidth + ELSE + grid((i - 1) * dimensions + j) = 0.0d0 + END IF + END DO + END DO + + ! Print the grid + print *, "Grid values:" + do i = 1, chunkLength + do j = 1, dimensions + print "(A,I4,A,F6.2)", "grid(", (i - 1) * dimensions + j, ") = ", grid((i - 1) * dimensions + j) + end do end do - ! finalize - call precicef_finalize() - write(*,*) "Exit FluidSolver" -end program FluidSolver + CALL precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) + + ! Check if Initial Data is Required and Write if Necessary + CALL precicef_requires_initial_data(bool) + IF (bool == 1) THEN + WRITE (*,*) 'Fluid: Writing initial data' + CALL precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure) + END IF + + ! Initialize Simulation Time + t = 0.0d0 + WRITE (*,*) "Initialize preCICE..." + CALL precicef_initialize() + + ! Read Initial Cross-Section Length + CALL precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, 0.0d0, crossSectionLength) + + ! Copy Current Cross-Section Length to Old Array + crossSectionLength_old = crossSectionLength + + ! initialize such that mass conservation is fulfilled + DO i = 1, chunkLength + velocity_old(i) = vel_in_0 * crossSectionLength_old(1) / crossSectionLength_old(i) + END DO + + ! Initialize Output Counter + out_counter = 0 + + ! Print all arrays with 2 decimal places + print *, "Index | Pressure | Pressure_Old | CrossSection | CrossSection_Old | Velocity | Velocity_Old" + do i = 1, chunkLength + print "(I5, 3X, F8.2, 3X, F13.2, 3X, F13.2, 3X, F16.2, 3X, F8.2, 3X, F13.2)", & + i - 1, pressure(i), pressure_old(i), crossSectionLength(i), & + crossSectionLength_old(i), velocity(i), velocity_old(i) + end do + + ! Main Coupling Loop + CALL precicef_is_coupling_ongoing(ongoing) + DO WHILE (ongoing /= 0) + ! Check if Writing a Checkpoint is Required + CALL precicef_requires_writing_checkpoint(bool) + IF (bool == 1) THEN + WRITE (*,*) 'Fluid: Writing iteration checkpoint' + END IF + + ! Get Maximum Time Step Size from preCICE + CALL precicef_get_max_time_step_size(dt) + + ! Compute Fluid Solution + CALL fluidComputeSolutionSerial( & + velocity_old, pressure_old, crossSectionLength_old, & + crossSectionLength, & + t + dt, & ! used for inlet velocity + domainSize, & + kappa, & + dt, & ! tau + velocity, pressure, & ! resulting velocity pressure + info) + + + CALL precicef_advance(dt) + + CALL precicef_get_max_time_step_size(dt) + + CALL precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, dt, crossSectionLength) + + CALL precicef_requires_reading_checkpoint(bool) + IF (bool == 1) THEN + WRITE (*,*) 'Fluid: Reading iteration checkpoint' + ELSE + t = t + dt + + CALL write_vtk(t, out_counter, outputFilePrefix, chunkLength, grid, velocity, pressure, crossSectionLength) + crossSectionLength_old = crossSectionLength + pressure_old = pressure + velocity_old = velocity + + out_counter = out_counter + 1 + END IF + + ! Check if Coupling is Still Ongoing + CALL precicef_is_coupling_ongoing(ongoing) + END DO + + ! Finalize preCICE Interface + CALL precicef_finalize() + WRITE (*,*) 'Exiting FluidSolver' + + ! Deallocate Dynamically Allocated Arrays + DEALLOCATE(pressure) + DEALLOCATE(pressure_old) + DEALLOCATE(crossSectionLength) + DEALLOCATE(crossSectionLength_old) + DEALLOCATE(velocity) + DEALLOCATE(velocity_old) + DEALLOCATE(grid) + DEALLOCATE(vertexIDs) + +END PROGRAM FluidSolver diff --git a/elastic-tube-1d/solid-fortran/CMakeLists.txt b/elastic-tube-1d/solid-fortran/CMakeLists.txt index ba56c7c3d..f063b730a 100644 --- a/elastic-tube-1d/solid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/solid-fortran/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION 3.10) -project(ElasticTube_Solid_Fortran LANGUAGES C Fortran DESCRIPTION "preCICE Fortran Solid Solver") +project(ElasticTube LANGUAGES Fortran C) if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE) @@ -9,13 +9,11 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) endif() message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}") - -find_package(precice REQUIRED CONFIG) - +find_package(precice 3 REQUIRED CONFIG) add_executable(SolidSolver + src/SolidComputeSolution.f90 src/SolidSolver.f90 ) - target_link_libraries(SolidSolver PRIVATE precice::precice) diff --git a/elastic-tube-1d/solid-fortran/clean.sh b/elastic-tube-1d/solid-fortran/clean.sh index bec0a3f61..ae7a54924 100755 --- a/elastic-tube-1d/solid-fortran/clean.sh +++ b/elastic-tube-1d/solid-fortran/clean.sh @@ -1,4 +1,7 @@ #!/usr/bin/env bash set -e -u -rm -rf build || true +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . diff --git a/elastic-tube-1d/solid-fortran/run.sh b/elastic-tube-1d/solid-fortran/run.sh index e362cf904..5263549b9 100755 --- a/elastic-tube-1d/solid-fortran/run.sh +++ b/elastic-tube-1d/solid-fortran/run.sh @@ -1,6 +1,9 @@ #!/usr/bin/env bash set -e -u +. ../../tools/log.sh || true +exec > >(tee --append "$LOGFILE") 2>&1 || true + if [ ! -d build ]; then mkdir build cd build @@ -10,3 +13,5 @@ if [ ! -d build ]; then fi ./build/SolidSolver ../precice-config.xml + +close_log diff --git a/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 b/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 index e69de29bb..a9540b6dc 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 @@ -0,0 +1,28 @@ +module SolidComputeSolution_mod + implicit none +contains + + subroutine SolidComputeSolution(chunkLength, pressure, crossSectionLength) + integer, intent(in) :: chunkLength + real(8), intent(in) :: pressure(1:chunkLength) + real(8), intent(inout) :: crossSectionLength(1:chunkLength) + + real(8) :: PI, E, r0, c_mk, c_mk2 + real(8) :: pressure0 + integer :: i + + PI = 3.14159265359d0 + E = 10000.d0 + r0 = 1.d0 / sqrt(PI) + c_mk = sqrt(E / (2.d0*r0)) + c_mk2 = c_mk * c_mk + pressure0 = 0.d0 + + do i=1, chunkLength + crossSectionLength(i) = ( (pressure0 - 2.d0 * c_mk2)**2 ) / & + ( (pressure(i) - 2.d0 * c_mk2)**2 ) + end do + + end subroutine SolidComputeSolution + +end module SolidComputeSolution_mod diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 8f4d4f23f..9db03ae13 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -1,42 +1,117 @@ -program SolidSolver - implicit none +PROGRAM SolidSolver + USE SolidComputeSolution_mod, ONLY: SolidComputeSolution + IMPLICIT NONE - character(len=512) :: configFileName - character(len=*), parameter :: solverName = 'Solid' + INTEGER, PARAMETER :: DP = KIND(1.0D0) - integer :: rank, commsize, ongoing - double precision :: dt + CHARACTER(LEN=512) :: configFileName + CHARACTER(LEN=256) :: solverName + CHARACTER(LEN=256) :: meshName, crossSectionLengthName, pressureName + INTEGER :: rank, commsize, ongoing, dimensions, bool + INTEGER :: domainSize, chunkLength + INTEGER :: i, j + INTEGER, ALLOCATABLE :: vertexIDs(:) + DOUBLE PRECISION, ALLOCATABLE :: pressure(:), crossSectionLength(:) + DOUBLE PRECISION, ALLOCATABLE :: grid(:) + DOUBLE PRECISION :: dt, tubeLength, dx - external :: precicef_create - external :: precicef_initialize - external :: precicef_is_coupling_ongoing - external :: precicef_get_max_time_step_size - external :: precicef_advance - external :: precicef_finalize + WRITE (*,*) 'Starting Solid Solver...' - write(*,*) "Starting minimal Solid Solver..." + IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN + WRITE (*,*) 'Solid: Usage: SolidSolver ' + WRITE (*,*) '' + STOP -1 + END IF - if(command_argument_count() /= 1) then - write(*,*) "Usage: SolidSolver " - stop - end if - call get_command_argument(1, configFileName) + CALL get_command_argument(1, configFileName) - rank = 0 + domainSize = 4 + chunkLength = domainSize + 1 + tubeLength = 10.0D0 + + WRITE (*,*) 'N: ', domainSize + WRITE (*,*) 'inputs: ', COMMAND_ARGUMENT_COUNT() + + solverName = 'Solid' + meshName = 'Solid-Nodes-Mesh' + crossSectionLengthName = 'CrossSectionLength' + pressureName = 'Pressure' + + rank = 0 commsize = 1 + CALL precicef_create(solverName, configFileName, rank, commsize) + WRITE (*,*) 'preCICE configured...' + + CALL precicef_get_mesh_dimensions(meshName, dimensions) + + ! Allocate Arrays + ALLOCATE(pressure(chunkLength)) + ALLOCATE(crossSectionLength(chunkLength)) + ALLOCATE(grid(dimensions * chunkLength)) + ALLOCATE(vertexIDs(chunkLength)) - call precicef_create(solverName, configFileName, rank, commsize) - - call precicef_initialize() + pressure = 0.0D0 + crossSectionLength = 1.0D0 + dx = tubeLength / domainSize + DO i = 1, chunkLength + grid((i - 1) * dimensions + 1) = dx * REAL(i - 1, DP) ! x-coordinate + grid((i - 1) * dimensions + 2) = 0.0D0 ! y-coordinate + vertexIDs(i) = i - 1 ! 0-based indexing here + END DO - call precicef_is_coupling_ongoing(ongoing) - do while(ongoing /= 0) - call precicef_get_max_time_step_size(dt) - call precicef_advance(dt) - call precicef_is_coupling_ongoing(ongoing) + ! Print the grid + print *, "Grid values:" + do i = 1, chunkLength + do j = 1, dimensions + print "(A,I4,A,F6.2)", "grid(", (i - 1) * dimensions + j, ") = ", grid((i - 1) * dimensions + j) + end do end do - call precicef_finalize() + CALL precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) + + ! Check if Initial Data is Required and Write if Necessary + CALL precicef_requires_initial_data(bool) + IF (bool == 1) THEN + CALL precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) + END IF + + WRITE (*,*) 'Initialize preCICE...' + CALL precicef_initialize() + + ! Coupling Loop + CALL precicef_is_coupling_ongoing(ongoing) + DO WHILE (ongoing /= 0) + + CALL precicef_requires_writing_checkpoint(bool) + IF (bool == 1) THEN + WRITE (*,*) 'Solid: Writing iteration checkpoint (not implemented).' + END IF + + CALL precicef_get_max_time_step_size(dt) + + CALL precicef_read_data(meshName, pressureName, chunkLength, vertexIDs, dt, pressure) + + CALL SolidComputeSolution(chunkLength, pressure, crossSectionLength) + + CALL precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) + + CALL precicef_advance(dt) + + CALL precicef_requires_reading_checkpoint(bool) + IF (bool == 1) THEN + WRITE (*,*) 'Solid: Reading iteration checkpoint (not implemented).' + END IF + + CALL precicef_is_coupling_ongoing(ongoing) + END DO + + WRITE (*,*) 'Exiting SolidSolver' + + CALL precicef_finalize() + + DEALLOCATE(pressure) + DEALLOCATE(crossSectionLength) + DEALLOCATE(grid) + DEALLOCATE(vertexIDs) - write(*,*) "Exit SolidSolver" -end program SolidSolver +END PROGRAM SolidSolver From 3fb4efcfd9452f2ff2ff0479c4adbc8836e3a308 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 17:04:01 +0100 Subject: [PATCH 04/24] Add output folder --- elastic-tube-1d/fluid-fortran/output/.keepme | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 elastic-tube-1d/fluid-fortran/output/.keepme diff --git a/elastic-tube-1d/fluid-fortran/output/.keepme b/elastic-tube-1d/fluid-fortran/output/.keepme new file mode 100644 index 000000000..e69de29bb From 05646918931c8f01c9494959cfbeaa20e7d1a662 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 19:02:42 +0100 Subject: [PATCH 05/24] Implement vtk writing, Add missing write_data call in fluid solver --- .../src/FluidComputeSolution.f90 | 12 +-- .../fluid-fortran/src/FluidSolver.f90 | 16 +++- .../fluid-fortran/src/utilities.f90 | 88 +++++++++++++++++++ .../solid-fortran/src/SolidSolver.f90 | 4 +- 4 files changed, 109 insertions(+), 11 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 index c2fd73996..f7aab2304 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 @@ -113,6 +113,7 @@ subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & if ((norm < tolerance .and. k > 1) .or. k > max_iterations) then exit + print *, "exiting" end if ! Initialize the LHS matrix @@ -180,14 +181,15 @@ subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & ! do i = 1, size(LHS, 1) ! print *, "LHS(", i, ",", i, ") =", LHS(i, i) ! end do - print *, "LHS Matrix (size: ", size(LHS, 1), "x", size(LHS, 2), "):" - do i = 1, size(LHS, 1) - write(*, '(10F8.2)') (LHS(i, j), j = 1, size(LHS, 2)) - end do + ! print *, "LHS Matrix (size: ", size(LHS, 1), "x", size(LHS, 2), "):" + ! do i = 1, size(LHS, 1) + ! write(*, '(10F8.2)') (LHS(i, j), j = 1, size(LHS, 2)) + ! end do call dgesv(nlhs, nrhs, LHS, nlhs, ipiv, Res, nlhs, info) if (info /= 0) then - return + ! return + print *, "Couldnt solve" end if ! Update velocity and pressure diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index e8c8c2768..633667bb4 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -1,6 +1,6 @@ PROGRAM FluidSolver USE FluidComputeSolution, ONLY: fluidComputeSolutionSerial - USE Utilities_mod, ONLY: write_vtk + USE utilities, ONLY: write_vtk IMPLICIT NONE ! Variable Declarations @@ -116,7 +116,6 @@ PROGRAM FluidSolver CALL precicef_requires_initial_data(bool) IF (bool == 1) THEN WRITE (*,*) 'Fluid: Writing initial data' - CALL precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure) END IF ! Initialize Simulation Time @@ -151,7 +150,7 @@ PROGRAM FluidSolver DO WHILE (ongoing /= 0) ! Check if Writing a Checkpoint is Required CALL precicef_requires_writing_checkpoint(bool) - IF (bool == 1) THEN + IF (bool.EQ.1) THEN WRITE (*,*) 'Fluid: Writing iteration checkpoint' END IF @@ -169,7 +168,16 @@ PROGRAM FluidSolver velocity, pressure, & ! resulting velocity pressure info) + ! Print all arrays with 2 decimal places + print *, "Index | Pressure | Pressure_Old | CrossSection | CrossSection_Old | Velocity | Velocity_Old" + do i = 1, chunkLength + print "(I5, 3X, F8.2, 3X, F13.2, 3X, F13.2, 3X, F16.2, 3X, F8.2, 3X, F13.2)", & + i - 1, pressure(i), pressure_old(i), crossSectionLength(i), & + crossSectionLength_old(i), velocity(i), velocity_old(i) + end do + CALL precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure) + CALL precicef_advance(dt) CALL precicef_get_max_time_step_size(dt) @@ -177,7 +185,7 @@ PROGRAM FluidSolver CALL precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, dt, crossSectionLength) CALL precicef_requires_reading_checkpoint(bool) - IF (bool == 1) THEN + IF (bool.EQ.1) THEN WRITE (*,*) 'Fluid: Reading iteration checkpoint' ELSE t = t + dt diff --git a/elastic-tube-1d/fluid-fortran/src/utilities.f90 b/elastic-tube-1d/fluid-fortran/src/utilities.f90 index e69de29bb..01fc489d6 100644 --- a/elastic-tube-1d/fluid-fortran/src/utilities.f90 +++ b/elastic-tube-1d/fluid-fortran/src/utilities.f90 @@ -0,0 +1,88 @@ +module utilities + implicit none + integer, parameter :: dp = kind(1.0d0) +contains + + subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & + grid, velocity, pressure, diameter) + implicit none + + ! Input Parameters + real(dp), intent(in) :: t + integer, intent(in) :: iteration + character(len=*), intent(in) :: filenamePrefix + integer, intent(in) :: nSlices + real(dp), intent(in) :: grid(:) + real(dp), intent(in) :: velocity(:) + real(dp), intent(in) :: pressure(:) + real(dp), intent(in) :: diameter(:) + + ! Local Variables + integer :: ioUnit, i, ioStatus + character(len=256) :: filename + + ! Generate Filename + write(filename, '(A,"_",I0,".vtk")') trim(filenamePrefix), iteration + + ! Print Status Message + print *, 'Writing timestep at t=', t, ' to ', trim(filename) + + ! Open File + open(newunit=ioUnit, file=trim(filename), status="replace", & + action="write", form="formatted", iostat=ioStatus) + if (ioStatus /= 0) then + print *, 'Error: Unable to open file ', trim(filename) + return + end if + + ! Write VTK Header + write(ioUnit, '(A)') '# vtk DataFile Version 2.0' + write(ioUnit, '(A)') '' + write(ioUnit, '(A)') 'ASCII' + write(ioUnit, '(A)') '' + write(ioUnit, '(A)') 'DATASET UNSTRUCTURED_GRID' + write(ioUnit, '(A)') '' + + ! Write POINTS Section + write(ioUnit, '(A,I0,A)') 'POINTS ', nSlices, ' float' + write(ioUnit, '(A)') '' + + ! Export Mesh Points + do i = 1, nSlices + write(ioUnit, '(F24.16,1X,F24.16,1X,F24.16)') grid(2*(i-1)+1), grid(2*(i-1)+2), 0.0_dp + end do + write(ioUnit, '(A)') '' + + ! Write POINT_DATA Section + write(ioUnit, '(A,I0)') 'POINT_DATA ', nSlices + write(ioUnit, '(A)') '' + + ! Export Velocity Vectors + write(ioUnit, '(A,A,A)') 'VECTORS ', 'velocity', ' float' + do i = 1, nSlices + write(ioUnit, '(F24.16,1X,F24.16,1X,F24.16)') velocity(i), 0.0_dp, 0.0_dp + end do + write(ioUnit, '(A)') '' + + ! Export Pressure Scalars + write(ioUnit, '(A,A,A)') 'SCALARS ', 'pressure', ' float' + write(ioUnit, '(A)') 'LOOKUP_TABLE default' + do i = 1, nSlices + write(ioUnit, '(F24.16)') pressure(i) + end do + write(ioUnit, '(A)') '' + + ! Export Diameter Scalars + write(ioUnit, '(A,A,A)') 'SCALARS ', 'diameter', ' float' + write(ioUnit, '(A)') 'LOOKUP_TABLE default' + do i = 1, nSlices + write(ioUnit, '(F24.16)') diameter(i) + end do + write(ioUnit, '(A)') '' + + ! Close File + close(ioUnit) + + end subroutine write_vtk + +end module utilities diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 9db03ae13..30d118e23 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -83,7 +83,7 @@ PROGRAM SolidSolver DO WHILE (ongoing /= 0) CALL precicef_requires_writing_checkpoint(bool) - IF (bool == 1) THEN + IF (bool.EQ.1) THEN WRITE (*,*) 'Solid: Writing iteration checkpoint (not implemented).' END IF @@ -98,7 +98,7 @@ PROGRAM SolidSolver CALL precicef_advance(dt) CALL precicef_requires_reading_checkpoint(bool) - IF (bool == 1) THEN + IF (bool.EQ.1) THEN WRITE (*,*) 'Solid: Reading iteration checkpoint (not implemented).' END IF From 05a60fa047668f7cc2b35e584bdf53f420e68a52 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 19:10:23 +0100 Subject: [PATCH 06/24] Set domainSize to 100 --- elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 | 2 +- elastic-tube-1d/solid-fortran/src/SolidSolver.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index 633667bb4..b7f519e71 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -48,7 +48,7 @@ PROGRAM FluidSolver pressureName = "Pressure" crossSectionLengthName = "CrossSectionLength" - domainSize = 4 + domainSize = 100 chunkLength = domainSize + 1 kappa = 100.0d0 L = 10.0d0 diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 30d118e23..2a81d40e6 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -25,7 +25,7 @@ PROGRAM SolidSolver CALL get_command_argument(1, configFileName) - domainSize = 4 + domainSize = 100 chunkLength = domainSize + 1 tubeLength = 10.0D0 From 3b6f0879a9cd42adebe2a49f80e9b5c5383c4a2b Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 19:38:28 +0100 Subject: [PATCH 07/24] Use ES format to output in scientific notation --- .../fluid-fortran/src/utilities.f90 | 153 ++++++++---------- 1 file changed, 67 insertions(+), 86 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/utilities.f90 b/elastic-tube-1d/fluid-fortran/src/utilities.f90 index 01fc489d6..15ca1fe9f 100644 --- a/elastic-tube-1d/fluid-fortran/src/utilities.f90 +++ b/elastic-tube-1d/fluid-fortran/src/utilities.f90 @@ -1,88 +1,69 @@ -module utilities - implicit none - integer, parameter :: dp = kind(1.0d0) -contains +MODULE Utilities + IMPLICIT NONE + INTEGER, PARAMETER :: dp = KIND(1.0D0) +CONTAINS - subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & + SUBROUTINE write_vtk(t, iteration, filenamePrefix, nSlices, & grid, velocity, pressure, diameter) - implicit none - - ! Input Parameters - real(dp), intent(in) :: t - integer, intent(in) :: iteration - character(len=*), intent(in) :: filenamePrefix - integer, intent(in) :: nSlices - real(dp), intent(in) :: grid(:) - real(dp), intent(in) :: velocity(:) - real(dp), intent(in) :: pressure(:) - real(dp), intent(in) :: diameter(:) - - ! Local Variables - integer :: ioUnit, i, ioStatus - character(len=256) :: filename - - ! Generate Filename - write(filename, '(A,"_",I0,".vtk")') trim(filenamePrefix), iteration - - ! Print Status Message - print *, 'Writing timestep at t=', t, ' to ', trim(filename) - - ! Open File - open(newunit=ioUnit, file=trim(filename), status="replace", & - action="write", form="formatted", iostat=ioStatus) - if (ioStatus /= 0) then - print *, 'Error: Unable to open file ', trim(filename) - return - end if - - ! Write VTK Header - write(ioUnit, '(A)') '# vtk DataFile Version 2.0' - write(ioUnit, '(A)') '' - write(ioUnit, '(A)') 'ASCII' - write(ioUnit, '(A)') '' - write(ioUnit, '(A)') 'DATASET UNSTRUCTURED_GRID' - write(ioUnit, '(A)') '' - - ! Write POINTS Section - write(ioUnit, '(A,I0,A)') 'POINTS ', nSlices, ' float' - write(ioUnit, '(A)') '' - - ! Export Mesh Points - do i = 1, nSlices - write(ioUnit, '(F24.16,1X,F24.16,1X,F24.16)') grid(2*(i-1)+1), grid(2*(i-1)+2), 0.0_dp - end do - write(ioUnit, '(A)') '' - - ! Write POINT_DATA Section - write(ioUnit, '(A,I0)') 'POINT_DATA ', nSlices - write(ioUnit, '(A)') '' - - ! Export Velocity Vectors - write(ioUnit, '(A,A,A)') 'VECTORS ', 'velocity', ' float' - do i = 1, nSlices - write(ioUnit, '(F24.16,1X,F24.16,1X,F24.16)') velocity(i), 0.0_dp, 0.0_dp - end do - write(ioUnit, '(A)') '' - - ! Export Pressure Scalars - write(ioUnit, '(A,A,A)') 'SCALARS ', 'pressure', ' float' - write(ioUnit, '(A)') 'LOOKUP_TABLE default' - do i = 1, nSlices - write(ioUnit, '(F24.16)') pressure(i) - end do - write(ioUnit, '(A)') '' - - ! Export Diameter Scalars - write(ioUnit, '(A,A,A)') 'SCALARS ', 'diameter', ' float' - write(ioUnit, '(A)') 'LOOKUP_TABLE default' - do i = 1, nSlices - write(ioUnit, '(F24.16)') diameter(i) - end do - write(ioUnit, '(A)') '' - - ! Close File - close(ioUnit) - - end subroutine write_vtk - -end module utilities + IMPLICIT NONE + + DOUBLE PRECISION, INTENT(IN) :: t + INTEGER, INTENT(IN) :: iteration + CHARACTER(LEN=*), INTENT(IN) :: filenamePrefix + INTEGER, INTENT(IN) :: nSlices + DOUBLE PRECISION, INTENT(IN) :: grid(:), velocity(:), pressure(:), diameter(:) + + INTEGER :: ioUnit, i, ioStatus + CHARACTER(LEN=256) :: filename + + WRITE(filename, '(A,"_",I0,".vtk")') TRIM(filenamePrefix), iteration + PRINT *, 'Writing timestep at t=', t, ' to ', TRIM(filename) + + OPEN(newunit=ioUnit, file=TRIM(filename), status="replace", action="write", form="formatted", iostat=ioStatus) + IF (ioStatus /= 0) THEN + PRINT *, 'Error: Unable to open file ', TRIM(filename) + RETURN + END IF + + WRITE(ioUnit, '(A)') '# vtk DataFile Version 2.0' + WRITE(ioUnit, '(A)') '' + WRITE(ioUnit, '(A)') 'ASCII' + WRITE(ioUnit, '(A)') '' + WRITE(ioUnit, '(A)') 'DATASET UNSTRUCTURED_GRID' + WRITE(ioUnit, '(A)') '' + + WRITE(ioUnit, '(A,I0,A)') 'POINTS ', nSlices, ' float' + WRITE(ioUnit, '(A)') '' + DO i = 1, nSlices + WRITE(ioUnit, '(ES24.16,1X,ES24.16,1X,ES24.16)') grid(2*(i-1)+1), grid(2*(i-1)+2), 0.0D0 + END DO + WRITE(ioUnit, '(A)') '' + + WRITE(ioUnit, '(A,I0)') 'POINT_DATA ', nSlices + WRITE(ioUnit, '(A)') '' + + WRITE(ioUnit, '(A,A,A)') 'VECTORS ', 'velocity', ' float' + DO i = 1, nSlices + WRITE(ioUnit, '(ES24.16,1X,ES24.16,1X,ES24.16)') velocity(i), 0.0D0, 0.0D0 + END DO + WRITE(ioUnit, '(A)') '' + + WRITE(ioUnit, '(A,A,A)') 'SCALARS ', 'pressure', ' float' + WRITE(ioUnit, '(A)') 'LOOKUP_TABLE default' + DO i = 1, nSlices + WRITE(ioUnit, '(ES24.16)') pressure(i) + END DO + WRITE(ioUnit, '(A)') '' + + WRITE(ioUnit, '(A,A,A)') 'SCALARS ', 'diameter', ' float' + WRITE(ioUnit, '(A)') 'LOOKUP_TABLE default' + DO i = 1, nSlices + WRITE(ioUnit, '(ES24.16)') diameter(i) + END DO + WRITE(ioUnit, '(A)') '' + + CLOSE(ioUnit) + + END SUBROUTINE write_vtk + +END MODULE Utilities From 1ad6396a7d601b9a106d71e43e42ebdb386d1bd1 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 21:58:31 +0100 Subject: [PATCH 08/24] Refactor: clean code, use snake_case for subroutines, apply dp for consistent precision --- .../src/FluidComputeSolution.f90 | 215 ++++++------- .../fluid-fortran/src/FluidSolver.f90 | 292 +++++++++--------- .../fluid-fortran/src/utilities.f90 | 108 +++---- .../src/SolidComputeSolution.f90 | 36 ++- .../solid-fortran/src/SolidSolver.f90 | 176 +++++------ 5 files changed, 403 insertions(+), 424 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 index f7aab2304..709f00bac 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 @@ -1,49 +1,46 @@ module FluidComputeSolution - use, intrinsic :: iso_c_binding, only: c_double, c_int implicit none - private - public :: fluidComputeSolutionSerial - + integer, parameter :: dp = kind(1.0d0) + contains - subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & - crossSectionLength_old, crossSectionLength, t, N, kappa, tau, & + subroutine fluid_compute_solution(velocity_old, pressure_old, & + crossSectionLength_old, crossSectionLength, t, n, kappa, tau, & velocity, pressure, info) - implicit none - - ! Function arguments - real(8), intent(in) :: velocity_old(:), pressure_old(:) - real(8), intent(in) :: crossSectionLength_old(:), crossSectionLength(:) - real(8), intent(in) :: t - integer, intent(in) :: N - real(8), intent(in) :: kappa, tau - real(8), intent(inout) :: velocity(:), pressure(:) + + real(dp), intent(in) :: velocity_old(:), pressure_old(:) + real(dp), intent(in) :: crossSectionLength_old(:), crossSectionLength(:) + real(dp), intent(in) :: t + integer, intent(in) :: n + real(dp), intent(in) :: kappa, tau + real(dp), intent(inout) :: velocity(:), pressure(:) integer, intent(out) :: info ! Local variables integer :: i, k, j - real(8), parameter :: PI = 3.141592653589793 - real(8), parameter :: E = 10000.0 - real(8), parameter :: c_mk2 = E / 2.0 * sqrt(PI) - real(8), parameter :: u0 = 10.0, ampl = 3.0, frequency = 10.0, & - t_shift = 0.0 - real(8), parameter :: tolerance = 1.0e-15 + real(dp), parameter :: pi = 3.141592653589793_dp + real(dp), parameter :: e = 10000.0_dp + real(dp), parameter :: c_mk2 = e / 2.0_dp * sqrt(pi) + real(dp), parameter :: u0 = 10.0_dp, ampl = 3.0_dp, frequency = 10.0_dp, & + t_shift = 0.0_dp + real(dp), parameter :: tolerance = 1.0e-15_dp integer, parameter :: max_iterations = 50 - real(8) :: alpha, L, dx, velocity_in, tmp2, norm_1, norm_2, norm - + real(dp) :: alpha, l, dx, velocity_in, tmp2, norm_1, norm_2, norm + ! LAPACK Variables integer :: nlhs, nrhs - real(8), allocatable :: Res(:) - real(8), allocatable :: LHS(:,:) + real(dp), allocatable :: res(:) + real(dp), allocatable :: lhs(:, :) integer, allocatable :: ipiv(:) - nlhs = 2 * N + 2 + nlhs = 2*N + 2 nrhs = 1 + ! Allocate arrays - allocate(Res(2 * N + 2)) - allocate(LHS(2 * N + 2, 2 * N + 2)) - allocate(ipiv(nlhs)) + allocate (Res(2*N + 2)) + allocate (LHS(2*N + 2, 2*N + 2)) + allocate (ipiv(nlhs)) velocity = velocity_old pressure = pressure_old @@ -51,13 +48,10 @@ subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & ! Stabilization intensity alpha = 0.0 !(N * kappa * tau) / (N * tau + 1); L = 10.0 - dx = L / kappa !1.0 / (N * kappa); + dx = L/kappa !1.0 / (N * kappa); - ! Initialize result variable + ! result variable info = 0 - print *, "CrossSectionLength:" - write(*, '(10F12.6)') crossSectionLength - print *, "Alpha:", alpha ! Nonlinear solver loop do k = 1, max_iterations @@ -67,49 +61,49 @@ subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & ! Compute residuals do i = 2, N ! Adjusted for 1-based indexing ! Momentum - Res(i) = (velocity_old(i) * crossSectionLength_old(i) - velocity(i) * crossSectionLength(i)) * dx / tau - Res(i) = Res(i) + 0.25 * (-crossSectionLength(i+1) * velocity(i) * velocity(i+1) - & - crossSectionLength(i) * velocity(i) * velocity(i+1)) - Res(i) = Res(i) + 0.25 * (-crossSectionLength(i+1) * velocity(i)**2 - & - crossSectionLength(i) * velocity(i)**2 + & - crossSectionLength(i) * velocity(i-1) * velocity(i) + & - crossSectionLength(i-1) * velocity(i-1) * velocity(i)) - Res(i) = Res(i) + 0.25 * (crossSectionLength(i-1) * velocity(i-1)**2 + & - crossSectionLength(i) * velocity(i-1)**2) - Res(i) = Res(i) + 0.25 * (crossSectionLength(i-1) * pressure(i-1) + & - crossSectionLength(i) * pressure(i-1) - & - crossSectionLength(i-1) * pressure(i) + & - crossSectionLength(i+1) * pressure(i) - & - crossSectionLength(i) * pressure(i+1) - & - crossSectionLength(i+1) * pressure(i+1)) + Res(i) = (velocity_old(i)*crossSectionLength_old(i) - velocity(i)*crossSectionLength(i))*dx/tau + Res(i) = Res(i) + 0.25*(-crossSectionLength(i + 1)*velocity(i)*velocity(i + 1) - & + crossSectionLength(i)*velocity(i)*velocity(i + 1)) + Res(i) = Res(i) + 0.25*(-crossSectionLength(i + 1)*velocity(i)**2 - & + crossSectionLength(i)*velocity(i)**2 + & + crossSectionLength(i)*velocity(i - 1)*velocity(i) + & + crossSectionLength(i - 1)*velocity(i - 1)*velocity(i)) + Res(i) = Res(i) + 0.25*(crossSectionLength(i - 1)*velocity(i - 1)**2 + & + crossSectionLength(i)*velocity(i - 1)**2) + Res(i) = Res(i) + 0.25*(crossSectionLength(i - 1)*pressure(i - 1) + & + crossSectionLength(i)*pressure(i - 1) - & + crossSectionLength(i - 1)*pressure(i) + & + crossSectionLength(i + 1)*pressure(i) - & + crossSectionLength(i)*pressure(i + 1) - & + crossSectionLength(i + 1)*pressure(i + 1)) ! Continuity - Res(i+N+1) = (crossSectionLength_old(i) - crossSectionLength(i)) * dx / tau - Res(i+N+1) = Res(i+N+1) + 0.25 * (crossSectionLength(i-1) * velocity(i-1) + & - crossSectionLength(i) * velocity(i-1) + & - crossSectionLength(i-1) * velocity(i) - & - crossSectionLength(i+1) * velocity(i) - & - crossSectionLength(i) * velocity(i+1) - & - crossSectionLength(i+1) * velocity(i+1)) - Res(i+N+1) = Res(i+N+1) + alpha * (pressure(i-1) - 2.0 * pressure(i) + pressure(i+1)) + Res(i + N + 1) = (crossSectionLength_old(i) - crossSectionLength(i))*dx/tau + Res(i + N + 1) = Res(i + N + 1) + 0.25*(crossSectionLength(i - 1)*velocity(i - 1) + & + crossSectionLength(i)*velocity(i - 1) + & + crossSectionLength(i - 1)*velocity(i) - & + crossSectionLength(i + 1)*velocity(i) - & + crossSectionLength(i)*velocity(i + 1) - & + crossSectionLength(i + 1)*velocity(i + 1)) + Res(i + N + 1) = Res(i + N + 1) + alpha*(pressure(i - 1) - 2.0*pressure(i) + pressure(i + 1)) end do ! Boundary conditions - velocity_in = u0 + ampl * sin(frequency * (t + t_shift) * PI) + velocity_in = u0 + ampl*sin(frequency*(t + t_shift)*PI) Res(1) = velocity_in - velocity(1) ! Pressure Inlet is linearly interpolated - Res(N+2) = -pressure(1) + 2.0 * pressure(2) - pressure(3) + Res(N + 2) = -pressure(1) + 2.0*pressure(2) - pressure(3) ! Velocity Outlet is linearly interpolated - Res(N+1) = -velocity(N+1) + 2.0 * velocity(N) - velocity(N-1) + Res(N + 1) = -velocity(N + 1) + 2.0*velocity(N) - velocity(N - 1) ! Pressure Outlet is "non-reflecting" - tmp2 = sqrt(c_mk2 - pressure_old(N+1) / 2.0) - & - (velocity(N+1) - velocity_old(N+1)) / 4.0 - Res(2*N+2) = -pressure(N+1) + 2.0 * (c_mk2 - tmp2**2) + tmp2 = sqrt(c_mk2 - pressure_old(N + 1)/2.0) - & + (velocity(N + 1) - velocity_old(N + 1))/4.0 + Res(2*N + 2) = -pressure(N + 1) + 2.0*(c_mk2 - tmp2**2) ! Compute residual norm norm_1 = sqrt(sum(Res**2)) norm_2 = sqrt(sum(pressure**2) + sum(velocity**2)) - norm = norm_1 / norm_2 + norm = norm_1/norm_2 if ((norm < tolerance .and. k > 1) .or. k > max_iterations) then exit @@ -122,69 +116,53 @@ subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & ! Populate LHS matrix do i = 2, N ! Momentum, Velocity - LHS(i, i-1) = LHS(i, i-1) + 0.25 * (-2.0 * crossSectionLength(i-1) * velocity(i-1) - & - 2.0 * crossSectionLength(i) * velocity(i-1) - & - crossSectionLength(i) * velocity(i) - crossSectionLength(i-1) * velocity(i)) - LHS(i, i) = LHS(i, i) + crossSectionLength(i) * dx / tau + & - 0.25 * (crossSectionLength(i+1) * velocity(i+1) + & - crossSectionLength(i) * velocity(i+1) + & - 2.0 * crossSectionLength(i+1) * velocity(i) + & - 2.0 * crossSectionLength(i) * velocity(i) - & - crossSectionLength(i) * velocity(i-1) - crossSectionLength(i-1) * velocity(i-1)) - LHS(i, i+1) = LHS(i, i+1) + 0.25 * (crossSectionLength(i+1) * velocity(i) + & - crossSectionLength(i) * velocity(i)) + LHS(i, i - 1) = LHS(i, i - 1) + 0.25*(-2.0*crossSectionLength(i - 1)*velocity(i - 1) - & + 2.0*crossSectionLength(i)*velocity(i - 1) - & + crossSectionLength(i)*velocity(i) - crossSectionLength(i - 1)*velocity(i)) + LHS(i, i) = LHS(i, i) + crossSectionLength(i)*dx/tau + & + 0.25*(crossSectionLength(i + 1)*velocity(i + 1) + & + crossSectionLength(i)*velocity(i + 1) + & + 2.0*crossSectionLength(i + 1)*velocity(i) + & + 2.0*crossSectionLength(i)*velocity(i) - & + crossSectionLength(i)*velocity(i - 1) - crossSectionLength(i - 1)*velocity(i - 1)) + LHS(i, i + 1) = LHS(i, i + 1) + 0.25*(crossSectionLength(i + 1)*velocity(i) + & + crossSectionLength(i)*velocity(i)) ! Momentum, Pressure - LHS(i, N+1+i-1) = LHS(i, N+1+i-1) - 0.25 * crossSectionLength(i-1) - & - 0.25 * crossSectionLength(i) - LHS(i, N+1+i) = LHS(i, N+1+i) + 0.25 * crossSectionLength(i-1) - & - 0.25 * crossSectionLength(i+1) - LHS(i, N+1+i+1) = LHS(i, N+1+i+1) + 0.25 * crossSectionLength(i) + & - 0.25 * crossSectionLength(i+1) + LHS(i, N + 1 + i - 1) = LHS(i, N + 1 + i - 1) - 0.25*crossSectionLength(i - 1) - & + 0.25*crossSectionLength(i) + LHS(i, N + 1 + i) = LHS(i, N + 1 + i) + 0.25*crossSectionLength(i - 1) - & + 0.25*crossSectionLength(i + 1) + LHS(i, N + 1 + i + 1) = LHS(i, N + 1 + i + 1) + 0.25*crossSectionLength(i) + & + 0.25*crossSectionLength(i + 1) ! Continuity, Velocity - LHS(i+N+1, i-1) = LHS(i+N+1, i-1) - 0.25 * crossSectionLength(i-1) - & - 0.25 * crossSectionLength(i) - LHS(i+N+1, i) = LHS(i+N+1, i) - 0.25 * crossSectionLength(i-1) + & - 0.25 * crossSectionLength(i+1) - LHS(i+N+1, i+1) = LHS(i+N+1, i+1) + 0.25 * crossSectionLength(i) + & - 0.25 * crossSectionLength(i+1) + LHS(i + N + 1, i - 1) = LHS(i + N + 1, i - 1) - 0.25*crossSectionLength(i - 1) - & + 0.25*crossSectionLength(i) + LHS(i + N + 1, i) = LHS(i + N + 1, i) - 0.25*crossSectionLength(i - 1) + & + 0.25*crossSectionLength(i + 1) + LHS(i + N + 1, i + 1) = LHS(i + N + 1, i + 1) + 0.25*crossSectionLength(i) + & + 0.25*crossSectionLength(i + 1) ! Continuity, Pressure - LHS(i+N+1, N+1+i-1) = LHS(i+N+1, N+1+i-1) - alpha - LHS(i+N+1, N+1+i) = LHS(i+N+1, N+1+i) + 2.0 * alpha - LHS(i+N+1, N+1+i+1) = LHS(i+N+1, N+1+i+1) - alpha + LHS(i + N + 1, N + 1 + i - 1) = LHS(i + N + 1, N + 1 + i - 1) - alpha + LHS(i + N + 1, N + 1 + i) = LHS(i + N + 1, N + 1 + i) + 2.0*alpha + LHS(i + N + 1, N + 1 + i + 1) = LHS(i + N + 1, N + 1 + i + 1) - alpha end do - ! Boundary conditions in LHS ! Velocity Inlet is prescribed LHS(1, 1) = 1.0 ! Pressure Inlet is linearly interpolated - LHS(N+2, N+2) = 1.0 - LHS(N+2, N+3) = -2.0 - LHS(N+2, N+4) = 1.0 + LHS(N + 2, N + 2) = 1.0 + LHS(N + 2, N + 3) = -2.0 + LHS(N + 2, N + 4) = 1.0 ! Velocity Outlet is linearly interpolated - LHS(N+1, N+1) = 1.0 - LHS(N+1, N) = -2.0 - LHS(N+1, N-1) = 1.0 + LHS(N + 1, N + 1) = 1.0 + LHS(N + 1, N) = -2.0 + LHS(N + 1, N - 1) = 1.0 ! Pressure Outlet is Non-Reflecting - LHS(2*N+2, 2*N+2) = 1.0 - LHS(2*N+2, N+1) = -(sqrt(c_mk2 - pressure_old(N+1) / 2.0) - (velocity(N+1) - velocity_old(N+1)) / 4.0) - - - ! Solve Ax = b using LAPACK - ! print *, "LHS Matrix (size: ", size(LHS, 1), "x", size(LHS, 2), "):" - ! do i = 1, size(LHS, 1) - ! write(*, '(10F12.6)') (LHS(i, j), j = 1, size(LHS, 2)) - ! end do - ! print *, "Diagonal elements of LHS:" - ! do i = 1, size(LHS, 1) - ! print *, "LHS(", i, ",", i, ") =", LHS(i, i) - ! end do - ! print *, "LHS Matrix (size: ", size(LHS, 1), "x", size(LHS, 2), "):" - ! do i = 1, size(LHS, 1) - ! write(*, '(10F8.2)') (LHS(i, j), j = 1, size(LHS, 2)) - ! end do + LHS(2*N + 2, 2*N + 2) = 1.0 + LHS(2*N + 2, N + 1) = -(sqrt(c_mk2 - pressure_old(N + 1)/2.0) - (velocity(N + 1) - velocity_old(N + 1))/4.0) call dgesv(nlhs, nrhs, LHS, nlhs, ipiv, Res, nlhs, info) if (info /= 0) then @@ -193,13 +171,14 @@ subroutine fluidComputeSolutionSerial(velocity_old, pressure_old, & end if ! Update velocity and pressure - do i = 1, N+1 + do i = 1, N + 1 velocity(i) = velocity(i) + Res(i) - pressure(i) = pressure(i) + Res(i+N+1) + pressure(i) = pressure(i) + Res(i + N + 1) end do end do ! Deallocate arrays deallocate(Res, LHS, ipiv) - end subroutine fluidComputeSolutionSerial -end module FluidComputeSolution \ No newline at end of file + + end subroutine fluid_compute_solution +end module FluidComputeSolution diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index b7f519e71..749311328 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -1,87 +1,82 @@ -PROGRAM FluidSolver - USE FluidComputeSolution, ONLY: fluidComputeSolutionSerial - USE utilities, ONLY: write_vtk - IMPLICIT NONE - - ! Variable Declarations - CHARACTER(LEN=512) :: configFileName - CHARACTER(LEN=50) :: solverName - CHARACTER(LEN=50) :: meshName, pressureName, crossSectionLengthName - CHARACTER(LEN=256) :: outputFilePrefix - INTEGER :: rank, commsize, ongoing, dimensions, bool - INTEGER :: domainSize, chunkLength - INTEGER :: i, j, info - DOUBLE PRECISION :: dt, t, cellwidth - DOUBLE PRECISION, ALLOCATABLE :: pressure(:), pressure_old(:) - DOUBLE PRECISION, ALLOCATABLE :: crossSectionLength(:), crossSectionLength_old(:) - DOUBLE PRECISION, ALLOCATABLE :: velocity(:), velocity_old(:) - INTEGER, ALLOCATABLE :: vertexIDs(:) - INTEGER :: out_counter - DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793d0 - DOUBLE PRECISION :: kappa, L - DOUBLE PRECISION :: r0, a0, u0, ampl, frequency, t_shift, p0, vel_in_0 - DOUBLE PRECISION, ALLOCATABLE :: grid(:) - - ! Start of Program - WRITE (*,*) 'Fluid: Starting Fortran solver...' - - ! Command-Line Argument Parsing - IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN - WRITE (*,*) "" - WRITE (*,*) "Fluid: Usage: FluidSolver " - STOP -1 - END IF - - CALL getarg(1, configFileName) +program FluidSolver + use FluidComputeSolution, only: fluid_compute_solution + use Utilities, only: write_vtk + implicit none + integer, parameter :: dp = kind(1.0d0) ! Double precision + + ! Variable declarations + character(LEN=512) :: configFileName + character(LEN=50) :: solverName + character(LEN=50) :: meshName, pressureName, crossSectionLengthName + character(LEN=256) :: outputFilePrefix + integer :: rank, commsize, ongoing, dimensions, bool + integer :: domainSize, chunkLength + integer :: i, j, info + real(dp) :: dt, t, cellwidth + real(dp), allocatable :: pressure(:), pressure_old(:) + real(dp), allocatable :: crossSectionLength(:), crossSectionLength_old(:) + real(dp), allocatable :: velocity(:), velocity_old(:) + integer, allocatable :: vertexIDs(:) + integer :: out_counter + real(dp), parameter :: pi = 3.141592653589793_dp + real(dp) :: kappa, l + real(dp) :: r0, a0, u0, ampl, frequency, t_shift, p0, vel_in_0 + real(dp), allocatable :: grid(:) + + + write(*, *) 'Fluid: Starting Fortran solver...' + + if (command_argument_count() /= 1) then + write(*, *) "" + write(*, *) "Fluid: Usage: FluidSolver " + stop -1 + end if + + call getarg(1, configFileName) solverName = 'Fluid' outputFilePrefix = './output/out_fluid' - ! Initialize preCICE Interface + ! Configure precice rank = 0 commsize = 1 - CALL precicef_create(solverName, configFileName, rank, commsize) - WRITE (*,*) "preCICE configured..." + call precicef_create(solverName, configFileName, rank, commsize) + write(*, *) "preCICE configured..." - ! Define Mesh and Data Names + ! Define mesh and data names meshName = "Fluid-Nodes-Mesh" pressureName = "Pressure" crossSectionLengthName = "CrossSectionLength" domainSize = 100 chunkLength = domainSize + 1 - kappa = 100.0d0 - L = 10.0d0 - - ! Get Mesh Dimensions from preCICE - CALL precicef_get_mesh_dimensions(meshName, dimensions) - - ! Allocate Arrays - ALLOCATE(vertexIDs(chunkLength)) - ALLOCATE(pressure(chunkLength)) - ALLOCATE(pressure_old(chunkLength)) - ALLOCATE(crossSectionLength(chunkLength)) - ALLOCATE(crossSectionLength_old(chunkLength)) - ALLOCATE(velocity(chunkLength)) - ALLOCATE(velocity_old(chunkLength)) - ALLOCATE(grid(dimensions * chunkLength)) - - ! Initialize vertexIDs (0-based IDs) - DO i = 1, chunkLength - vertexIDs(i) = i - 1 - END DO - - ! Initialize Physical Parameters - r0 = 1.0d0 / SQRT(PI) - a0 = r0**2 * PI - u0 = 10.0d0 - ampl = 3.0d0 - frequency = 10.0d0 - t_shift = 0.0d0 - p0 = 0.0d0 - vel_in_0 = u0 + ampl * SIN(frequency * (t_shift) * PI) - - ! Initialize Data Arrays + kappa = 100.0_dp + l = 10.0_dp + + ! Get mesh dimensions + call precicef_get_mesh_dimensions(meshName, dimensions) + + ! Allocate arrays + allocate(vertexIDs(chunkLength)) + allocate(pressure(chunkLength)) + allocate(pressure_old(chunkLength)) + allocate(crossSectionLength(chunkLength)) + allocate(crossSectionLength_old(chunkLength)) + allocate(velocity(chunkLength)) + allocate(velocity_old(chunkLength)) + allocate(grid(dimensions*chunkLength)) + + ! Initialize physical parameters + r0 = 1.0_dp / sqrt(pi) + a0 = r0**2 * pi + u0 = 10.0_dp + ampl = 3.0_dp + frequency = 10.0_dp + t_shift = 0.0_dp + p0 = 0.0_dp + vel_in_0 = u0 + ampl * sin(frequency * (t_shift) * pi) + + ! Initialize data arrays pressure = p0 pressure_old = pressure crossSectionLength = a0 @@ -89,52 +84,55 @@ PROGRAM FluidSolver velocity = vel_in_0 velocity_old = velocity - ! Initialize Grid Coordinates - cellwidth = L / REAL(domainSize, KIND=8) - DO i = 1, chunkLength - DO j = 1, dimensions - IF (j == 1) THEN - grid((i - 1) * dimensions + j) = REAL(i - 1, KIND=8) * cellwidth - ELSE - grid((i - 1) * dimensions + j) = 0.0d0 - END IF - END DO - END DO + ! Initialize grid coordinates + cellwidth = l / real(domainSize, dp) + do i = 1, chunkLength + do j = 1, dimensions + if (j == 1) then + grid((i - 1)*dimensions + j) = real(i - 1, dp) * cellwidth + else + grid((i - 1)*dimensions + j) = 0.0_dp + end if + end do + end do + + ! Initialize vertexIDs (0-based IDs) + do i = 1, chunkLength + vertexIDs(i) = i - 1 + end do ! Print the grid print *, "Grid values:" do i = 1, chunkLength do j = 1, dimensions - print "(A,I4,A,F6.2)", "grid(", (i - 1) * dimensions + j, ") = ", grid((i - 1) * dimensions + j) + print "(A,I4,A,F6.2)", "grid(", (i - 1)*dimensions + j, ") = ", grid((i - 1)*dimensions + j) end do end do - - CALL precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) + call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) ! Check if Initial Data is Required and Write if Necessary - CALL precicef_requires_initial_data(bool) - IF (bool == 1) THEN - WRITE (*,*) 'Fluid: Writing initial data' - END IF + call precicef_requires_initial_data(bool) + if (bool == 1) then + write (*, *) 'Fluid: Writing initial data' + end if - ! Initialize Simulation Time - t = 0.0d0 - WRITE (*,*) "Initialize preCICE..." - CALL precicef_initialize() + + t = 0.0d0 + write (*, *) "Initialize preCICE..." + call precicef_initialize() - ! Read Initial Cross-Section Length - CALL precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, 0.0d0, crossSectionLength) + ! read initial cross-Section length + call precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, 0.0d0, crossSectionLength) - ! Copy Current Cross-Section Length to Old Array + ! Copy current cross-Section length to old array crossSectionLength_old = crossSectionLength ! initialize such that mass conservation is fulfilled - DO i = 1, chunkLength - velocity_old(i) = vel_in_0 * crossSectionLength_old(1) / crossSectionLength_old(i) - END DO + do i = 1, chunkLength + velocity_old(i) = vel_in_0*crossSectionLength_old(1)/crossSectionLength_old(i) + end do - ! Initialize Output Counter out_counter = 0 ! Print all arrays with 2 decimal places @@ -143,30 +141,29 @@ PROGRAM FluidSolver print "(I5, 3X, F8.2, 3X, F13.2, 3X, F13.2, 3X, F16.2, 3X, F8.2, 3X, F13.2)", & i - 1, pressure(i), pressure_old(i), crossSectionLength(i), & crossSectionLength_old(i), velocity(i), velocity_old(i) - end do - - ! Main Coupling Loop - CALL precicef_is_coupling_ongoing(ongoing) - DO WHILE (ongoing /= 0) - ! Check if Writing a Checkpoint is Required - CALL precicef_requires_writing_checkpoint(bool) - IF (bool.EQ.1) THEN - WRITE (*,*) 'Fluid: Writing iteration checkpoint' - END IF - - ! Get Maximum Time Step Size from preCICE - CALL precicef_get_max_time_step_size(dt) - - ! Compute Fluid Solution - CALL fluidComputeSolutionSerial( & + end do + + ! Main coupling loop + call precicef_is_coupling_ongoing(ongoing) + do while (ongoing /= 0) + ! Check if writing a checkpoint is required + call precicef_requires_writing_checkpoint(bool) + if (bool .eq. 1) then + write (*, *) 'Fluid: Writing iteration checkpoint' + end if + + call precicef_get_max_time_step_size(dt) + + ! solve + call fluid_compute_solution( & velocity_old, pressure_old, crossSectionLength_old, & crossSectionLength, & t + dt, & ! used for inlet velocity - domainSize, & - kappa, & + domainSize, & + kappa, & dt, & ! tau velocity, pressure, & ! resulting velocity pressure - info) + info) ! Print all arrays with 2 decimal places print *, "Index | Pressure | Pressure_Old | CrossSection | CrossSection_Old | Velocity | Velocity_Old" @@ -174,46 +171,45 @@ PROGRAM FluidSolver print "(I5, 3X, F8.2, 3X, F13.2, 3X, F13.2, 3X, F16.2, 3X, F8.2, 3X, F13.2)", & i - 1, pressure(i), pressure_old(i), crossSectionLength(i), & crossSectionLength_old(i), velocity(i), velocity_old(i) - end do + end do - CALL precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure) + call precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure) - CALL precicef_advance(dt) + call precicef_advance(dt) - CALL precicef_get_max_time_step_size(dt) + call precicef_get_max_time_step_size(dt) - CALL precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, dt, crossSectionLength) + call precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, dt, crossSectionLength) - CALL precicef_requires_reading_checkpoint(bool) - IF (bool.EQ.1) THEN - WRITE (*,*) 'Fluid: Reading iteration checkpoint' - ELSE + call precicef_requires_reading_checkpoint(bool) + if (bool .eq. 1) then + write (*, *) 'Fluid: Reading iteration checkpoint' + else t = t + dt - CALL write_vtk(t, out_counter, outputFilePrefix, chunkLength, grid, velocity, pressure, crossSectionLength) + call write_vtk(t, out_counter, outputFilePrefix, chunkLength, grid, velocity, pressure, crossSectionLength) crossSectionLength_old = crossSectionLength pressure_old = pressure velocity_old = velocity out_counter = out_counter + 1 - END IF + end if ! Check if Coupling is Still Ongoing - CALL precicef_is_coupling_ongoing(ongoing) - END DO - - ! Finalize preCICE Interface - CALL precicef_finalize() - WRITE (*,*) 'Exiting FluidSolver' - - ! Deallocate Dynamically Allocated Arrays - DEALLOCATE(pressure) - DEALLOCATE(pressure_old) - DEALLOCATE(crossSectionLength) - DEALLOCATE(crossSectionLength_old) - DEALLOCATE(velocity) - DEALLOCATE(velocity_old) - DEALLOCATE(grid) - DEALLOCATE(vertexIDs) - -END PROGRAM FluidSolver + call precicef_is_coupling_ongoing(ongoing) + end do + + ! finalize precice and deallocate arrays + call precicef_finalize() + write (*, *) 'Exiting FluidSolver' + + deallocate(pressure) + deallocate(pressure_old) + deallocate(crossSectionLength) + deallocate(crossSectionLength_old) + deallocate(velocity) + deallocate(velocity_old) + deallocate(grid) + deallocate(vertexIDs) + +end program FluidSolver diff --git a/elastic-tube-1d/fluid-fortran/src/utilities.f90 b/elastic-tube-1d/fluid-fortran/src/utilities.f90 index 15ca1fe9f..4cf04561b 100644 --- a/elastic-tube-1d/fluid-fortran/src/utilities.f90 +++ b/elastic-tube-1d/fluid-fortran/src/utilities.f90 @@ -1,69 +1,69 @@ -MODULE Utilities - IMPLICIT NONE - INTEGER, PARAMETER :: dp = KIND(1.0D0) -CONTAINS +module Utilities + implicit none + integer, parameter :: dp = kind(1.0D0) +contains - SUBROUTINE write_vtk(t, iteration, filenamePrefix, nSlices, & + subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & grid, velocity, pressure, diameter) - IMPLICIT NONE + implicit none - DOUBLE PRECISION, INTENT(IN) :: t - INTEGER, INTENT(IN) :: iteration - CHARACTER(LEN=*), INTENT(IN) :: filenamePrefix - INTEGER, INTENT(IN) :: nSlices - DOUBLE PRECISION, INTENT(IN) :: grid(:), velocity(:), pressure(:), diameter(:) + double precision, intent(IN) :: t + integer, intent(IN) :: iteration + character(LEN=*), intent(IN) :: filenamePrefix + integer, intent(IN) :: nSlices + double precision, intent(IN) :: grid(:), velocity(:), pressure(:), diameter(:) - INTEGER :: ioUnit, i, ioStatus - CHARACTER(LEN=256) :: filename + integer :: ioUnit, i, ioStatus + character(LEN=256) :: filename - WRITE(filename, '(A,"_",I0,".vtk")') TRIM(filenamePrefix), iteration - PRINT *, 'Writing timestep at t=', t, ' to ', TRIM(filename) + write (filename, '(A,"_",I0,".vtk")') trim(filenamePrefix), iteration + print *, 'Writing timestep at t=', t, ' to ', trim(filename) - OPEN(newunit=ioUnit, file=TRIM(filename), status="replace", action="write", form="formatted", iostat=ioStatus) - IF (ioStatus /= 0) THEN - PRINT *, 'Error: Unable to open file ', TRIM(filename) - RETURN - END IF + open (newunit=ioUnit, file=trim(filename), status="replace", action="write", form="formatted", iostat=ioStatus) + if (ioStatus /= 0) then + print *, 'Error: Unable to open file ', trim(filename) + return + end if - WRITE(ioUnit, '(A)') '# vtk DataFile Version 2.0' - WRITE(ioUnit, '(A)') '' - WRITE(ioUnit, '(A)') 'ASCII' - WRITE(ioUnit, '(A)') '' - WRITE(ioUnit, '(A)') 'DATASET UNSTRUCTURED_GRID' - WRITE(ioUnit, '(A)') '' + write (ioUnit, '(A)') '# vtk DataFile Version 2.0' + write (ioUnit, '(A)') '' + write (ioUnit, '(A)') 'ASCII' + write (ioUnit, '(A)') '' + write (ioUnit, '(A)') 'DATASET UNSTRUCTURED_GRID' + write (ioUnit, '(A)') '' - WRITE(ioUnit, '(A,I0,A)') 'POINTS ', nSlices, ' float' - WRITE(ioUnit, '(A)') '' - DO i = 1, nSlices - WRITE(ioUnit, '(ES24.16,1X,ES24.16,1X,ES24.16)') grid(2*(i-1)+1), grid(2*(i-1)+2), 0.0D0 - END DO - WRITE(ioUnit, '(A)') '' + write (ioUnit, '(A,I0,A)') 'POINTS ', nSlices, ' float' + write (ioUnit, '(A)') '' + do i = 1, nSlices + write (ioUnit, '(ES24.16,1X,ES24.16,1X,ES24.16)') grid(2*(i - 1) + 1), grid(2*(i - 1) + 2), 0.0D0 + end do + write (ioUnit, '(A)') '' - WRITE(ioUnit, '(A,I0)') 'POINT_DATA ', nSlices - WRITE(ioUnit, '(A)') '' + write (ioUnit, '(A,I0)') 'POINT_DATA ', nSlices + write (ioUnit, '(A)') '' - WRITE(ioUnit, '(A,A,A)') 'VECTORS ', 'velocity', ' float' - DO i = 1, nSlices - WRITE(ioUnit, '(ES24.16,1X,ES24.16,1X,ES24.16)') velocity(i), 0.0D0, 0.0D0 - END DO - WRITE(ioUnit, '(A)') '' + write (ioUnit, '(A,A,A)') 'VECTORS ', 'velocity', ' float' + do i = 1, nSlices + write (ioUnit, '(ES24.16,1X,ES24.16,1X,ES24.16)') velocity(i), 0.0D0, 0.0D0 + end do + write (ioUnit, '(A)') '' - WRITE(ioUnit, '(A,A,A)') 'SCALARS ', 'pressure', ' float' - WRITE(ioUnit, '(A)') 'LOOKUP_TABLE default' - DO i = 1, nSlices - WRITE(ioUnit, '(ES24.16)') pressure(i) - END DO - WRITE(ioUnit, '(A)') '' + write (ioUnit, '(A,A,A)') 'SCALARS ', 'pressure', ' float' + write (ioUnit, '(A)') 'LOOKUP_TABLE default' + do i = 1, nSlices + write (ioUnit, '(ES24.16)') pressure(i) + end do + write (ioUnit, '(A)') '' - WRITE(ioUnit, '(A,A,A)') 'SCALARS ', 'diameter', ' float' - WRITE(ioUnit, '(A)') 'LOOKUP_TABLE default' - DO i = 1, nSlices - WRITE(ioUnit, '(ES24.16)') diameter(i) - END DO - WRITE(ioUnit, '(A)') '' + write (ioUnit, '(A,A,A)') 'SCALARS ', 'diameter', ' float' + write (ioUnit, '(A)') 'LOOKUP_TABLE default' + do i = 1, nSlices + write (ioUnit, '(ES24.16)') diameter(i) + end do + write (ioUnit, '(A)') '' - CLOSE(ioUnit) + close (ioUnit) - END SUBROUTINE write_vtk + end subroutine write_vtk -END MODULE Utilities +end module Utilities diff --git a/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 b/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 index a9540b6dc..28aeb72b4 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidComputeSolution.f90 @@ -1,28 +1,32 @@ -module SolidComputeSolution_mod +module SolidComputeSolution implicit none + integer, parameter :: dp = kind(1.0d0) + contains - subroutine SolidComputeSolution(chunkLength, pressure, crossSectionLength) + subroutine solid_compute_solution(chunkLength, pressure, crossSectionLength) integer, intent(in) :: chunkLength - real(8), intent(in) :: pressure(1:chunkLength) - real(8), intent(inout) :: crossSectionLength(1:chunkLength) + real(dp), intent(in) :: pressure(1:chunkLength) + real(dp), intent(inout) :: crossSectionLength(1:chunkLength) - real(8) :: PI, E, r0, c_mk, c_mk2 - real(8) :: pressure0 + real(dp) :: pi, e, r0, c_mk, c_mk2 + real(dp) :: pressure0 integer :: i - PI = 3.14159265359d0 - E = 10000.d0 - r0 = 1.d0 / sqrt(PI) - c_mk = sqrt(E / (2.d0*r0)) + ! constants + pi = 3.141592653589793_dp + e = 10000.0_dp + r0 = 1.0_dp / sqrt(pi) + c_mk = sqrt(e / (2.0_dp * r0)) c_mk2 = c_mk * c_mk - pressure0 = 0.d0 + pressure0 = 0.0_dp - do i=1, chunkLength - crossSectionLength(i) = ( (pressure0 - 2.d0 * c_mk2)**2 ) / & - ( (pressure(i) - 2.d0 * c_mk2)**2 ) + ! Update crossSectionLength based on pressure + do i = 1, chunkLength + crossSectionLength(i) = ((pressure0 - 2.0_dp * c_mk2)**2) / & + ((pressure(i) - 2.0_dp * c_mk2)**2) end do - end subroutine SolidComputeSolution + end subroutine solid_compute_solution -end module SolidComputeSolution_mod +end module SolidComputeSolution diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 2a81d40e6..67746f532 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -1,117 +1,117 @@ -PROGRAM SolidSolver - USE SolidComputeSolution_mod, ONLY: SolidComputeSolution - IMPLICIT NONE - - INTEGER, PARAMETER :: DP = KIND(1.0D0) - - CHARACTER(LEN=512) :: configFileName - CHARACTER(LEN=256) :: solverName - CHARACTER(LEN=256) :: meshName, crossSectionLengthName, pressureName - INTEGER :: rank, commsize, ongoing, dimensions, bool - INTEGER :: domainSize, chunkLength - INTEGER :: i, j - INTEGER, ALLOCATABLE :: vertexIDs(:) - DOUBLE PRECISION, ALLOCATABLE :: pressure(:), crossSectionLength(:) - DOUBLE PRECISION, ALLOCATABLE :: grid(:) - DOUBLE PRECISION :: dt, tubeLength, dx - - WRITE (*,*) 'Starting Solid Solver...' - - IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN - WRITE (*,*) 'Solid: Usage: SolidSolver ' - WRITE (*,*) '' - STOP -1 - END IF - - CALL get_command_argument(1, configFileName) - - domainSize = 100 +program SolidSolver + use SolidComputeSolution, only: solid_compute_solution + implicit none + + integer, parameter :: dp = kind(1.0d0) + + character(len=512) :: configFileName + character(len=256) :: solverName + character(len=256) :: meshName, crossSectionLengthName, pressureName + integer :: rank, commsize, ongoing, dimensions, bool + integer :: domainSize, chunkLength + integer :: i, j + integer, allocatable :: vertexIDs(:) + real(dp), allocatable :: pressure(:), crossSectionLength(:) + real(dp), allocatable :: grid(:) + real(dp) :: dt, tubeLength, dx + + write(*, *) 'Starting Solid Solver...' + + if (command_argument_count() /= 1) then + write(*, *) 'Solid: Usage: SolidSolver ' + write(*, *) '' + stop -1 + end if + + call get_command_argument(1, configFileName) + + domainSize = 100 chunkLength = domainSize + 1 - tubeLength = 10.0D0 + tubeLength = 10.0_dp - WRITE (*,*) 'N: ', domainSize - WRITE (*,*) 'inputs: ', COMMAND_ARGUMENT_COUNT() + write(*, *) 'N: ', domainSize + write(*, *) 'inputs: ', command_argument_count() - solverName = 'Solid' - meshName = 'Solid-Nodes-Mesh' - crossSectionLengthName = 'CrossSectionLength' - pressureName = 'Pressure' + solverName = 'Solid' + meshName = 'Solid-Nodes-Mesh' + crossSectionLengthName = 'CrossSectionLength' + pressureName = 'Pressure' - rank = 0 + rank = 0 commsize = 1 - CALL precicef_create(solverName, configFileName, rank, commsize) - WRITE (*,*) 'preCICE configured...' - - CALL precicef_get_mesh_dimensions(meshName, dimensions) - - ! Allocate Arrays - ALLOCATE(pressure(chunkLength)) - ALLOCATE(crossSectionLength(chunkLength)) - ALLOCATE(grid(dimensions * chunkLength)) - ALLOCATE(vertexIDs(chunkLength)) - - pressure = 0.0D0 - crossSectionLength = 1.0D0 - dx = tubeLength / domainSize - DO i = 1, chunkLength - grid((i - 1) * dimensions + 1) = dx * REAL(i - 1, DP) ! x-coordinate - grid((i - 1) * dimensions + 2) = 0.0D0 ! y-coordinate + call precicef_create(solverName, configFileName, rank, commsize) + write(*, *) 'preCICE configured...' + + call precicef_get_mesh_dimensions(meshName, dimensions) + + ! Allocate arrays + allocate(pressure(chunkLength)) + allocate(crossSectionLength(chunkLength)) + allocate(grid(dimensions*chunkLength)) + allocate(vertexIDs(chunkLength)) + + pressure = 0.0_dp + crossSectionLength = 1.0_dp + dx = tubeLength / real(domainSize, dp) + do i = 1, chunkLength + grid((i - 1)*dimensions + 1) = dx * real(i - 1, dp) ! x-coordinate + grid((i - 1)*dimensions + 2) = 0.0_dp ! y-coordinate vertexIDs(i) = i - 1 ! 0-based indexing here - END DO + end do ! Print the grid print *, "Grid values:" do i = 1, chunkLength do j = 1, dimensions - print "(A,I4,A,F6.2)", "grid(", (i - 1) * dimensions + j, ") = ", grid((i - 1) * dimensions + j) + print "(A,I4,A,F6.2)", "grid(", (i - 1)*dimensions + j, ") = ", grid((i - 1)*dimensions + j) end do end do - CALL precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) + call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) - ! Check if Initial Data is Required and Write if Necessary - CALL precicef_requires_initial_data(bool) - IF (bool == 1) THEN - CALL precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) - END IF + ! Check if initial data is required and write if necessary + call precicef_requires_initial_data(bool) + if (bool == 1) then + call precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) + end if - WRITE (*,*) 'Initialize preCICE...' - CALL precicef_initialize() + write (*, *) 'Initialize preCICE...' + call precicef_initialize() - ! Coupling Loop - CALL precicef_is_coupling_ongoing(ongoing) - DO WHILE (ongoing /= 0) + ! Coupling loop + call precicef_is_coupling_ongoing(ongoing) + do while (ongoing /= 0) - CALL precicef_requires_writing_checkpoint(bool) - IF (bool.EQ.1) THEN - WRITE (*,*) 'Solid: Writing iteration checkpoint (not implemented).' - END IF + call precicef_requires_writing_checkpoint(bool) + if (bool .eq. 1) then + write (*, *) 'Solid: Writing iteration checkpoint (not implemented).' + end if - CALL precicef_get_max_time_step_size(dt) + call precicef_get_max_time_step_size(dt) - CALL precicef_read_data(meshName, pressureName, chunkLength, vertexIDs, dt, pressure) + call precicef_read_data(meshName, pressureName, chunkLength, vertexIDs, dt, pressure) - CALL SolidComputeSolution(chunkLength, pressure, crossSectionLength) + call solid_compute_solution(chunkLength, pressure, crossSectionLength) - CALL precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) + call precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) - CALL precicef_advance(dt) + call precicef_advance(dt) - CALL precicef_requires_reading_checkpoint(bool) - IF (bool.EQ.1) THEN - WRITE (*,*) 'Solid: Reading iteration checkpoint (not implemented).' - END IF + call precicef_requires_reading_checkpoint(bool) + if (bool .eq. 1) then + write (*, *) 'Solid: Reading iteration checkpoint (not implemented).' + end if - CALL precicef_is_coupling_ongoing(ongoing) - END DO + call precicef_is_coupling_ongoing(ongoing) + end do - WRITE (*,*) 'Exiting SolidSolver' + write (*, *) 'Exiting SolidSolver' - CALL precicef_finalize() + call precicef_finalize() - DEALLOCATE(pressure) - DEALLOCATE(crossSectionLength) - DEALLOCATE(grid) - DEALLOCATE(vertexIDs) + deallocate(pressure) + deallocate(crossSectionLength) + deallocate(grid) + deallocate(vertexIDs) -END PROGRAM SolidSolver +end program SolidSolver From 9975e75672b4e776f2048771a2444a4216b72444 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 22:29:22 +0100 Subject: [PATCH 09/24] Remove unused var --- elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 index 709f00bac..7c0981270 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 @@ -1,7 +1,7 @@ module FluidComputeSolution implicit none integer, parameter :: dp = kind(1.0d0) - + contains subroutine fluid_compute_solution(velocity_old, pressure_old, & @@ -17,7 +17,7 @@ subroutine fluid_compute_solution(velocity_old, pressure_old, & integer, intent(out) :: info ! Local variables - integer :: i, k, j + integer :: i, k real(dp), parameter :: pi = 3.141592653589793_dp real(dp), parameter :: e = 10000.0_dp real(dp), parameter :: c_mk2 = e / 2.0_dp * sqrt(pi) From 0738a43818a1f41e3a738aaf1fdb7f6ccebc918f Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 22:30:03 +0100 Subject: [PATCH 10/24] apply dp for consistent precision --- elastic-tube-1d/fluid-fortran/src/utilities.f90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/utilities.f90 b/elastic-tube-1d/fluid-fortran/src/utilities.f90 index 4cf04561b..e2587ae39 100644 --- a/elastic-tube-1d/fluid-fortran/src/utilities.f90 +++ b/elastic-tube-1d/fluid-fortran/src/utilities.f90 @@ -7,11 +7,11 @@ subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & grid, velocity, pressure, diameter) implicit none - double precision, intent(IN) :: t - integer, intent(IN) :: iteration - character(LEN=*), intent(IN) :: filenamePrefix - integer, intent(IN) :: nSlices - double precision, intent(IN) :: grid(:), velocity(:), pressure(:), diameter(:) + real(dp), intent(IN) :: t + integer, intent(IN) :: iteration + character(LEN=*), intent(IN) :: filenamePrefix + integer, intent(IN) :: nSlices + real(dp), intent(IN) :: grid(:), velocity(:), pressure(:), diameter(:) integer :: ioUnit, i, ioStatus character(LEN=256) :: filename @@ -25,6 +25,7 @@ subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & return end if + ! Write vtk headers write (ioUnit, '(A)') '# vtk DataFile Version 2.0' write (ioUnit, '(A)') '' write (ioUnit, '(A)') 'ASCII' @@ -32,6 +33,7 @@ subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & write (ioUnit, '(A)') 'DATASET UNSTRUCTURED_GRID' write (ioUnit, '(A)') '' + ! Write points write (ioUnit, '(A,I0,A)') 'POINTS ', nSlices, ' float' write (ioUnit, '(A)') '' do i = 1, nSlices @@ -42,12 +44,14 @@ subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & write (ioUnit, '(A,I0)') 'POINT_DATA ', nSlices write (ioUnit, '(A)') '' + ! Write velocity vector field write (ioUnit, '(A,A,A)') 'VECTORS ', 'velocity', ' float' do i = 1, nSlices write (ioUnit, '(ES24.16,1X,ES24.16,1X,ES24.16)') velocity(i), 0.0D0, 0.0D0 end do write (ioUnit, '(A)') '' + ! Write pressure write (ioUnit, '(A,A,A)') 'SCALARS ', 'pressure', ' float' write (ioUnit, '(A)') 'LOOKUP_TABLE default' do i = 1, nSlices @@ -55,6 +59,7 @@ subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & end do write (ioUnit, '(A)') '' + ! Write diameter write (ioUnit, '(A,A,A)') 'SCALARS ', 'diameter', ' float' write (ioUnit, '(A)') 'LOOKUP_TABLE default' do i = 1, nSlices From 1f285bb668a8e109a77399e2d856e2e3f2013f03 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 22:50:08 +0100 Subject: [PATCH 11/24] Remove debug prints --- .../src/FluidComputeSolution.f90 | 4 +-- .../fluid-fortran/src/FluidSolver.f90 | 36 ++++--------------- .../solid-fortran/src/SolidSolver.f90 | 12 ++----- 3 files changed, 9 insertions(+), 43 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 index 7c0981270..df9604816 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidComputeSolution.f90 @@ -107,7 +107,6 @@ subroutine fluid_compute_solution(velocity_old, pressure_old, & if ((norm < tolerance .and. k > 1) .or. k > max_iterations) then exit - print *, "exiting" end if ! Initialize the LHS matrix @@ -166,8 +165,7 @@ subroutine fluid_compute_solution(velocity_old, pressure_old, & call dgesv(nlhs, nrhs, LHS, nlhs, ipiv, Res, nlhs, info) if (info /= 0) then - ! return - print *, "Couldnt solve" + write(*, *) "Linear Solver not converged!, Info: ", info end if ! Update velocity and pressure diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index 749311328..7576a932e 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -101,14 +101,6 @@ program FluidSolver vertexIDs(i) = i - 1 end do - ! Print the grid - print *, "Grid values:" - do i = 1, chunkLength - do j = 1, dimensions - print "(A,I4,A,F6.2)", "grid(", (i - 1)*dimensions + j, ") = ", grid((i - 1)*dimensions + j) - end do - end do - call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) ! Check if Initial Data is Required and Write if Necessary @@ -117,8 +109,6 @@ program FluidSolver write (*, *) 'Fluid: Writing initial data' end if - - t = 0.0d0 write (*, *) "Initialize preCICE..." call precicef_initialize() @@ -133,23 +123,16 @@ program FluidSolver velocity_old(i) = vel_in_0*crossSectionLength_old(1)/crossSectionLength_old(i) end do + t = 0.0d0 out_counter = 0 - ! Print all arrays with 2 decimal places - print *, "Index | Pressure | Pressure_Old | CrossSection | CrossSection_Old | Velocity | Velocity_Old" - do i = 1, chunkLength - print "(I5, 3X, F8.2, 3X, F13.2, 3X, F13.2, 3X, F16.2, 3X, F8.2, 3X, F13.2)", & - i - 1, pressure(i), pressure_old(i), crossSectionLength(i), & - crossSectionLength_old(i), velocity(i), velocity_old(i) - end do - ! Main coupling loop call precicef_is_coupling_ongoing(ongoing) do while (ongoing /= 0) - ! Check if writing a checkpoint is required + ! checkpointing is required in implicit coupling call precicef_requires_writing_checkpoint(bool) if (bool .eq. 1) then - write (*, *) 'Fluid: Writing iteration checkpoint' + ! nothing end if call precicef_get_max_time_step_size(dt) @@ -165,14 +148,6 @@ program FluidSolver velocity, pressure, & ! resulting velocity pressure info) - ! Print all arrays with 2 decimal places - print *, "Index | Pressure | Pressure_Old | CrossSection | CrossSection_Old | Velocity | Velocity_Old" - do i = 1, chunkLength - print "(I5, 3X, F8.2, 3X, F13.2, 3X, F13.2, 3X, F16.2, 3X, F8.2, 3X, F13.2)", & - i - 1, pressure(i), pressure_old(i), crossSectionLength(i), & - crossSectionLength_old(i), velocity(i), velocity_old(i) - end do - call precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure) call precicef_advance(dt) @@ -183,8 +158,9 @@ program FluidSolver call precicef_requires_reading_checkpoint(bool) if (bool .eq. 1) then - write (*, *) 'Fluid: Reading iteration checkpoint' + ! not yet converged else + ! converged, advance in time t = t + dt call write_vtk(t, out_counter, outputFilePrefix, chunkLength, grid, velocity, pressure, crossSectionLength) @@ -195,7 +171,7 @@ program FluidSolver out_counter = out_counter + 1 end if - ! Check if Coupling is Still Ongoing + ! Check if coupling is still ongoing call precicef_is_coupling_ongoing(ongoing) end do diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 67746f532..77854b5b5 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -59,14 +59,6 @@ program SolidSolver vertexIDs(i) = i - 1 ! 0-based indexing here end do - ! Print the grid - print *, "Grid values:" - do i = 1, chunkLength - do j = 1, dimensions - print "(A,I4,A,F6.2)", "grid(", (i - 1)*dimensions + j, ") = ", grid((i - 1)*dimensions + j) - end do - end do - call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) ! Check if initial data is required and write if necessary @@ -84,7 +76,7 @@ program SolidSolver call precicef_requires_writing_checkpoint(bool) if (bool .eq. 1) then - write (*, *) 'Solid: Writing iteration checkpoint (not implemented).' + ! Do nothing here end if call precicef_get_max_time_step_size(dt) @@ -99,7 +91,7 @@ program SolidSolver call precicef_requires_reading_checkpoint(bool) if (bool .eq. 1) then - write (*, *) 'Solid: Reading iteration checkpoint (not implemented).' + ! nothing end if call precicef_is_coupling_ongoing(ongoing) From 82cffe9e620d9a9c23acca50d90d8d8e58ad273c Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 22:55:41 +0100 Subject: [PATCH 12/24] Remove comment --- elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index 7576a932e..a173580a6 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -4,7 +4,6 @@ program FluidSolver implicit none integer, parameter :: dp = kind(1.0d0) ! Double precision - ! Variable declarations character(LEN=512) :: configFileName character(LEN=50) :: solverName character(LEN=50) :: meshName, pressureName, crossSectionLengthName From 29f1ace87d2e1507bc3c4fbbef21cdfc546f564f Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Mon, 6 Jan 2025 23:51:00 +0100 Subject: [PATCH 13/24] Include Fortran watchpoint in plot-all script --- elastic-tube-1d/plot-all.sh | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/elastic-tube-1d/plot-all.sh b/elastic-tube-1d/plot-all.sh index 8ba00ee03..26975482d 100755 --- a/elastic-tube-1d/plot-all.sh +++ b/elastic-tube-1d/plot-all.sh @@ -12,15 +12,16 @@ gnuplot -p << EOF plot \ "fluid-cpp/precice-Fluid-watchpoint-Middle.log" using 1:4 with linespoints title "cpp", \ "fluid-python/precice-Fluid-watchpoint-Middle.log" using 1:4 with linespoints lw 2 title "python", \ - "fluid-rust/precice-Fluid-watchpoint-Middle.log" using 1:4 with linespoints title "rust" + "fluid-rust/precice-Fluid-watchpoint-Middle.log" using 1:4 with linespoints title "rust", \ + "fluid-fortran/precice-Fluid-watchpoint-Middle.log" using 1:4 with linespoints title "fortran" set title 'Pressure of elastic-tube at x=5' set ylabel 'pressure' plot \ "fluid-cpp/precice-Fluid-watchpoint-Middle.log" using 1:5 with linespoints title "cpp", \ "fluid-python/precice-Fluid-watchpoint-Middle.log" using 1:5 with linespoints title "python", \ - "fluid-rust/precice-Fluid-watchpoint-Middle.log" using 1:5 with linespoints title "rust" - - set xlabel 'time [s]' + "fluid-rust/precice-Fluid-watchpoint-Middle.log" using 1:5 with linespoints title "rust", \ + "fluid-fortran/precice-Fluid-watchpoint-Middle.log" using 1:5 with linespoints title "fortran" + set xlabel 'time [s]' EOF From 99b484acd115512486c16099afba963b133a10a2 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 00:00:49 +0100 Subject: [PATCH 14/24] Improve print statement formatting --- elastic-tube-1d/fluid-fortran/src/utilities.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastic-tube-1d/fluid-fortran/src/utilities.f90 b/elastic-tube-1d/fluid-fortran/src/utilities.f90 index e2587ae39..e844e0a64 100644 --- a/elastic-tube-1d/fluid-fortran/src/utilities.f90 +++ b/elastic-tube-1d/fluid-fortran/src/utilities.f90 @@ -17,7 +17,7 @@ subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & character(LEN=256) :: filename write (filename, '(A,"_",I0,".vtk")') trim(filenamePrefix), iteration - print *, 'Writing timestep at t=', t, ' to ', trim(filename) + print '(A, F7.6, A, A)', 'writing timestep at t=', t, ' to ', trim(filename) open (newunit=ioUnit, file=trim(filename), status="replace", action="write", form="formatted", iostat=ioStatus) if (ioStatus /= 0) then From a0d8609eb3d6c3fe23d63bee085ef7c9e058ff52 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 13:15:11 +0100 Subject: [PATCH 15/24] Remove debug compiler flags and clean up --- elastic-tube-1d/fluid-fortran/CMakeLists.txt | 75 -------------------- 1 file changed, 75 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/CMakeLists.txt b/elastic-tube-1d/fluid-fortran/CMakeLists.txt index e14fecfc4..248b02f77 100644 --- a/elastic-tube-1d/fluid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/fluid-fortran/CMakeLists.txt @@ -1,7 +1,6 @@ cmake_minimum_required(VERSION 3.10) project(ElasticTube LANGUAGES Fortran C) -# Set default build type to Debug if not specified if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE) set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS @@ -9,89 +8,15 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) endif() message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}") -# Find required packages find_package(precice 3 REQUIRED CONFIG) find_package(LAPACK REQUIRED) -# Add executable add_executable(FluidSolver src/FluidComputeSolution.f90 src/utilities.f90 src/FluidSolver.f90 ) -# Link libraries target_link_libraries(FluidSolver PRIVATE precice::precice) target_link_libraries(FluidSolver PRIVATE ${LAPACK_LIBRARIES} ${LAPACK_LINKER_FLAGS}) -# ============================ -# Add Compiler Flags for Debugging -# ============================ - -# Define a function to add compiler flags based on compiler ID -function(add_fortran_debug_flags target) - if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") - # GNU Fortran (gfortran) specific flags - target_compile_options(${target} PRIVATE - $<$:-fbounds-check> - $<$:-fbacktrace> - $<$:-fcheck=all> - $<$:-g> - $<$:-Wall> - $<$:-Wextra> - ) - elseif(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") - # Intel Fortran (ifort) specific flags - target_compile_options(${target} PRIVATE - $<$:-check bounds> - $<$:-g> - $<$:-warn all> - ) - elseif(CMAKE_Fortran_COMPILER_ID MATCHES "PGI" OR - CMAKE_Fortran_COMPILER_ID MATCHES "NAG") - # PGI, NAG, and other Fortran compilers - target_compile_options(${target} PRIVATE - $<$:-C> - $<$:-g> - $<$:-Wall> - ) - else() - message(WARNING "Compiler ${CMAKE_Fortran_COMPILER_ID} not explicitly supported for debug flags.") - # Generic debug flags - target_compile_options(${target} PRIVATE - $<$:-g> - $<$:-Wall> - ) - endif() -endfunction() - -# Apply the debug flags to FluidSolver -add_fortran_debug_flags(FluidSolver) - -# ============================ -# Optionally, Add Linker Flags (if needed) -# ============================ - -# For example, to link LAPACK, you might need to specify additional flags -# However, typically, find_package(LAPACK) handles this -# If you have specific linker flags, you can add them as follows: -# target_link_libraries(FluidSolver PRIVATE ${LAPACK_LIBRARIES} ${LAPACK_LINKER_FLAGS}) - -# ============================ -# Enable Fortran Module Generation (Optional) -# ============================ - -# If you have modules that need to be visible to other targets, set the appropriate properties -# For example: -# set_target_properties(FluidSolver PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/modules) - -# ============================ -# Set CMake Policies (Optional) -# ============================ - -# To ensure compatibility with newer CMake versions -# cmake_policy(SET CMP0074 NEW) - -# ============================ -# End of CMakeLists.txt -# ============================ From 7e228639d99a7abf91208445d8d729d75a7d8030 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 13:43:25 +0100 Subject: [PATCH 16/24] Set Fortran standard and enable warnings in debug builds --- elastic-tube-1d/fluid-fortran/CMakeLists.txt | 7 +++++++ elastic-tube-1d/solid-fortran/CMakeLists.txt | 6 ++++++ 2 files changed, 13 insertions(+) diff --git a/elastic-tube-1d/fluid-fortran/CMakeLists.txt b/elastic-tube-1d/fluid-fortran/CMakeLists.txt index 248b02f77..1978e7d9a 100644 --- a/elastic-tube-1d/fluid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/fluid-fortran/CMakeLists.txt @@ -1,6 +1,8 @@ cmake_minimum_required(VERSION 3.10) project(ElasticTube LANGUAGES Fortran C) +set(CMAKE_Fortran_STANDARD 2003) + if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE) set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS @@ -8,6 +10,11 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) endif() message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}") +# Add bounds check and warnings in debug build +if(CMAKE_BUILD_TYPE STREQUAL "Debug" AND CMAKE_Fortran_COMPILER_ID MATCHES "GNU") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wall -Wextra -fbounds-check") +endif() + find_package(precice 3 REQUIRED CONFIG) find_package(LAPACK REQUIRED) diff --git a/elastic-tube-1d/solid-fortran/CMakeLists.txt b/elastic-tube-1d/solid-fortran/CMakeLists.txt index f063b730a..c33915351 100644 --- a/elastic-tube-1d/solid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/solid-fortran/CMakeLists.txt @@ -2,6 +2,8 @@ cmake_minimum_required(VERSION 3.10) project(ElasticTube LANGUAGES Fortran C) +set(CMAKE_Fortran_STANDARD 2003) + if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Choose the type of build." FORCE) set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS @@ -9,6 +11,10 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) endif() message(STATUS "Build configuration: ${CMAKE_BUILD_TYPE}") +if(CMAKE_BUILD_TYPE STREQUAL "Debug" AND CMAKE_Fortran_COMPILER_ID MATCHES "GNU") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wall -Wextra -fbounds-check") +endif() + find_package(precice 3 REQUIRED CONFIG) add_executable(SolidSolver From c8ccc2a6399fb7ff8570ff34764ca14251782624 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 13:43:37 +0100 Subject: [PATCH 17/24] Remove unused variable --- elastic-tube-1d/solid-fortran/src/SolidSolver.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 77854b5b5..2a7b523e6 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -9,7 +9,7 @@ program SolidSolver character(len=256) :: meshName, crossSectionLengthName, pressureName integer :: rank, commsize, ongoing, dimensions, bool integer :: domainSize, chunkLength - integer :: i, j + integer :: i integer, allocatable :: vertexIDs(:) real(dp), allocatable :: pressure(:), crossSectionLength(:) real(dp), allocatable :: grid(:) From e4d602373fd58664cb05f24235d560523fce56c7 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 14:27:13 +0100 Subject: [PATCH 18/24] Set character lengths to 50 --- elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 | 2 +- elastic-tube-1d/solid-fortran/src/SolidSolver.f90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index a173580a6..f59b30df4 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -4,7 +4,7 @@ program FluidSolver implicit none integer, parameter :: dp = kind(1.0d0) ! Double precision - character(LEN=512) :: configFileName + character(LEN=50) :: configFileName character(LEN=50) :: solverName character(LEN=50) :: meshName, pressureName, crossSectionLengthName character(LEN=256) :: outputFilePrefix diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 2a7b523e6..9de1b1d63 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -4,9 +4,9 @@ program SolidSolver integer, parameter :: dp = kind(1.0d0) - character(len=512) :: configFileName - character(len=256) :: solverName - character(len=256) :: meshName, crossSectionLengthName, pressureName + character(len=50) :: configFileName + character(len=50) :: solverName + character(len=50) :: meshName, crossSectionLengthName, pressureName integer :: rank, commsize, ongoing, dimensions, bool integer :: domainSize, chunkLength integer :: i From 993fcff68d740b0143e0d9ea6dc8309223f78e06 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 14:34:19 +0100 Subject: [PATCH 19/24] Add preCICE Fortran module file (precice.f90) --- elastic-tube-1d/solid-fortran/src/precice.f90 | 315 ++++++++++++++++++ 1 file changed, 315 insertions(+) create mode 100644 elastic-tube-1d/solid-fortran/src/precice.f90 diff --git a/elastic-tube-1d/solid-fortran/src/precice.f90 b/elastic-tube-1d/solid-fortran/src/precice.f90 new file mode 100644 index 000000000..13f880e2b --- /dev/null +++ b/elastic-tube-1d/solid-fortran/src/precice.f90 @@ -0,0 +1,315 @@ +module precice + use, intrinsic :: iso_c_binding + implicit none + + interface + + subroutine precicef_create(participantName, configFileName, & + & solverProcessIndex, solverProcessSize, & + & participantNameLength, configFileNameLength) & + & bind(c, name='precicef_create_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: participantName + character(kind=c_char), dimension(*) :: configFileName + integer(kind=c_int) :: solverProcessIndex + integer(kind=c_int) :: solverProcessSize + integer(kind=c_int), value :: participantNameLength + integer(kind=c_int), value :: configFileNameLength + end subroutine precicef_create + + subroutine precicef_initialize() & + & bind(c, name='precicef_initialize_') + + use, intrinsic :: iso_c_binding + end subroutine precicef_initialize + + subroutine precicef_advance(timestepLengthLimit) & + & bind(c, name='precicef_advance_') + + use, intrinsic :: iso_c_binding + real(kind=c_double) :: timestepLengthLimit + end subroutine precicef_advance + + subroutine precicef_finalize() & + & bind(c, name='precicef_finalize_') + + use, intrinsic :: iso_c_binding + end subroutine precicef_finalize + + subroutine precicef_requires_reading_checkpoint(isRequired) & + & bind(c, name='precicef_requires_reading_checkpoint_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isRequired + end subroutine precicef_requires_reading_checkpoint + + subroutine precicef_requires_writing_checkpoint(isRequired) & + & bind(c, name='precicef_requires_writing_checkpoint_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isRequired + end subroutine precicef_requires_writing_checkpoint + + subroutine precicef_get_mesh_dimensions(meshName, dimensions, & + & meshNameLength) & + & bind(c, name='precicef_get_mesh_dimensions_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: dimensions + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_get_mesh_dimensions + + subroutine precicef_get_data_dimensions(meshName, dataName, & + & dimensions, meshNameLength, dataNameLength) & + & bind(c, name='precicef_get_data_dimensions_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: dimensions + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_get_data_dimensions + + subroutine precicef_is_coupling_ongoing(isOngoing) & + & bind(c, name='precicef_is_coupling_ongoing_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isOngoing + end subroutine precicef_is_coupling_ongoing + + subroutine precicef_is_time_window_complete(isComplete) & + & bind(c, name='precicef_is_time_window_complete_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isComplete + end subroutine precicef_is_time_window_complete + + subroutine precicef_get_max_time_step_size(maxTimeStepSize) & + & bind(c, name='precicef_get_max_time_step_size_') + + use, intrinsic :: iso_c_binding + real(kind=c_double) :: maxTimeStepSize + end subroutine precicef_get_max_time_step_size + + subroutine precicef_requires_mesh_connectivity_for(meshName, required, meshNameLength) & + & bind(c, name='precicef_requires_mesh_connectivity_for_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: required + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_requires_mesh_connectivity_for + + subroutine precicef_set_vertex(meshName, position, vertexID, meshNameLength) & + & bind(c, name='precicef_set_vertex_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + real(kind=c_double) :: coordinates(3) + integer(kind=c_int) :: id + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_vertex + + subroutine precicef_get_mesh_vertex_size(meshName, meshSize, meshNameLength) & + & bind(c, name='precicef_get_mesh_vertex_size_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: meshSize + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_get_mesh_vertex_size + + subroutine precicef_set_vertices(meshName, size, coordinates, ids, meshNameLength) & + & bind(c, name='precicef_set_vertices_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + real(kind=c_double) :: coordinates(*) + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_vertices + + subroutine precicef_set_edge(meshName, firstVertexID, secondVertexID, & + & meshNameLength) & + & bind(c, name='precicef_set_edge_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstVertexID + integer(kind=c_int) :: secondVertexID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_edge + + subroutine precicef_set_mesh_edges(meshName, size, ids, meshNameLength) & + & bind(c, name='precicef_set_mesh_edges_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_edges + + subroutine precicef_set_triangle(meshName, firstEdgeID, secondEdgeID, & + & thirdEdgeID, meshNameLength) & + & bind(c, name='precicef_set_triangle_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstEdgeID + integer(kind=c_int) :: secondEdgeID + integer(kind=c_int) :: thirdEdgeID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_triangle + + subroutine precicef_set_quad(meshName, firstVertexID, secondVertexID, & + & thirdVertexID, fourthVertexID, & + & meshNameLength ) & + & bind(c, name='precicef_set_quad_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstVertexID + integer(kind=c_int) :: secondVertexID + integer(kind=c_int) :: thirdVertexID + integer(kind=c_int) :: fourthVertexID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_quad + + subroutine precicef_set_mesh_quads(meshName, size, ids, meshNameLength) & + & bind(c, name='precicef_set_mesh_quads_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_quads + + subroutine precicef_set_tetrahedron(meshName, firstVertexID, secondVertexID, & + & thirdVertexID, fourthVertexID, & + & meshNameLength) & + & bind(c, name='precicef_set_tetrahedron_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstVertexID + integer(kind=c_int) :: secondVertexID + integer(kind=c_int) :: thirdVertexID + integer(kind=c_int) :: fourthVertexID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_tetrahedron + + subroutine precicef_set_mesh_tetrahedra(meshName, size, ids, meshNameLength) & + & bind(c, name='precicef_set_mesh_tetrahedra_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_tetrahedra + + subroutine precicef_requires_initial_data(isRequired) & + & bind(c, name='precicef_requires_initial_data_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isRequired + end subroutine precicef_requires_initial_data + + subroutine precicef_write_data(meshName, dataName, size, ids, & + & values, meshNameLength, dataNameLength) & + & bind(c, name='precicef_write_data_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: values(*) + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_write_data + + subroutine precicef_read_data(meshName, dataName, size, ids, & + & relativeReadTime, values, meshNameLength, & + & dataNameLength) & + & bind(c, name='precicef_read_data_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: relativeReadTime + real(kind=c_double) :: values(*) + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_read_data + + subroutine precicef_set_mesh_access_region(meshName, boundingBox, & + & meshNameLength) & + & bind(c, name='precicef_set_mesh_access_region_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + real(kind=c_double) :: boundingBox(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_access_region + + subroutine precicef_get_mesh_vertex_ids_and_coordinates(meshName, size, & + & ids, coordinates, & + & meshNameLength) & + & bind(c, name='precicef_get_mesh_vertex_ids_and_coordinates_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: coordinates(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_get_mesh_vertex_ids_and_coordinates + + subroutine precicef_requires_gradient_data_for(meshName, dataName, & + & required, meshNameLength, & + & dataNameLength) & + & bind(c, name='precicef_requires_gradient_data_for_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: required + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_requires_gradient_data_for + + subroutine precicef_write_gradient_data(meshName, dataName, size, ids, & + & gradients, meshNameLength, & + & dataNameLength) & + & bind(c, name='precicef_write_gradient_data_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: gradients(*) + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_write_gradient_data + + subroutine precicef_get_version_information(versionInfo, lengthVersionInfo) & + & bind(c, name="precicef_get_version_information_") + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: versionInfo + integer(kind=c_int), value :: lengthVersionInfo + end subroutine precicef_get_version_information + + end interface + +end module precice From 76dac6b38ad5f119bd903066889b2d4e5bc6da08 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 14:36:34 +0100 Subject: [PATCH 20/24] Add preCICE Fortran module file (precice.f90) --- elastic-tube-1d/fluid-fortran/src/precice.f90 | 315 ++++++++++++++++++ 1 file changed, 315 insertions(+) create mode 100644 elastic-tube-1d/fluid-fortran/src/precice.f90 diff --git a/elastic-tube-1d/fluid-fortran/src/precice.f90 b/elastic-tube-1d/fluid-fortran/src/precice.f90 new file mode 100644 index 000000000..13f880e2b --- /dev/null +++ b/elastic-tube-1d/fluid-fortran/src/precice.f90 @@ -0,0 +1,315 @@ +module precice + use, intrinsic :: iso_c_binding + implicit none + + interface + + subroutine precicef_create(participantName, configFileName, & + & solverProcessIndex, solverProcessSize, & + & participantNameLength, configFileNameLength) & + & bind(c, name='precicef_create_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: participantName + character(kind=c_char), dimension(*) :: configFileName + integer(kind=c_int) :: solverProcessIndex + integer(kind=c_int) :: solverProcessSize + integer(kind=c_int), value :: participantNameLength + integer(kind=c_int), value :: configFileNameLength + end subroutine precicef_create + + subroutine precicef_initialize() & + & bind(c, name='precicef_initialize_') + + use, intrinsic :: iso_c_binding + end subroutine precicef_initialize + + subroutine precicef_advance(timestepLengthLimit) & + & bind(c, name='precicef_advance_') + + use, intrinsic :: iso_c_binding + real(kind=c_double) :: timestepLengthLimit + end subroutine precicef_advance + + subroutine precicef_finalize() & + & bind(c, name='precicef_finalize_') + + use, intrinsic :: iso_c_binding + end subroutine precicef_finalize + + subroutine precicef_requires_reading_checkpoint(isRequired) & + & bind(c, name='precicef_requires_reading_checkpoint_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isRequired + end subroutine precicef_requires_reading_checkpoint + + subroutine precicef_requires_writing_checkpoint(isRequired) & + & bind(c, name='precicef_requires_writing_checkpoint_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isRequired + end subroutine precicef_requires_writing_checkpoint + + subroutine precicef_get_mesh_dimensions(meshName, dimensions, & + & meshNameLength) & + & bind(c, name='precicef_get_mesh_dimensions_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: dimensions + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_get_mesh_dimensions + + subroutine precicef_get_data_dimensions(meshName, dataName, & + & dimensions, meshNameLength, dataNameLength) & + & bind(c, name='precicef_get_data_dimensions_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: dimensions + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_get_data_dimensions + + subroutine precicef_is_coupling_ongoing(isOngoing) & + & bind(c, name='precicef_is_coupling_ongoing_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isOngoing + end subroutine precicef_is_coupling_ongoing + + subroutine precicef_is_time_window_complete(isComplete) & + & bind(c, name='precicef_is_time_window_complete_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isComplete + end subroutine precicef_is_time_window_complete + + subroutine precicef_get_max_time_step_size(maxTimeStepSize) & + & bind(c, name='precicef_get_max_time_step_size_') + + use, intrinsic :: iso_c_binding + real(kind=c_double) :: maxTimeStepSize + end subroutine precicef_get_max_time_step_size + + subroutine precicef_requires_mesh_connectivity_for(meshName, required, meshNameLength) & + & bind(c, name='precicef_requires_mesh_connectivity_for_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: required + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_requires_mesh_connectivity_for + + subroutine precicef_set_vertex(meshName, position, vertexID, meshNameLength) & + & bind(c, name='precicef_set_vertex_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + real(kind=c_double) :: coordinates(3) + integer(kind=c_int) :: id + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_vertex + + subroutine precicef_get_mesh_vertex_size(meshName, meshSize, meshNameLength) & + & bind(c, name='precicef_get_mesh_vertex_size_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: meshSize + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_get_mesh_vertex_size + + subroutine precicef_set_vertices(meshName, size, coordinates, ids, meshNameLength) & + & bind(c, name='precicef_set_vertices_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + real(kind=c_double) :: coordinates(*) + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_vertices + + subroutine precicef_set_edge(meshName, firstVertexID, secondVertexID, & + & meshNameLength) & + & bind(c, name='precicef_set_edge_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstVertexID + integer(kind=c_int) :: secondVertexID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_edge + + subroutine precicef_set_mesh_edges(meshName, size, ids, meshNameLength) & + & bind(c, name='precicef_set_mesh_edges_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_edges + + subroutine precicef_set_triangle(meshName, firstEdgeID, secondEdgeID, & + & thirdEdgeID, meshNameLength) & + & bind(c, name='precicef_set_triangle_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstEdgeID + integer(kind=c_int) :: secondEdgeID + integer(kind=c_int) :: thirdEdgeID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_triangle + + subroutine precicef_set_quad(meshName, firstVertexID, secondVertexID, & + & thirdVertexID, fourthVertexID, & + & meshNameLength ) & + & bind(c, name='precicef_set_quad_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstVertexID + integer(kind=c_int) :: secondVertexID + integer(kind=c_int) :: thirdVertexID + integer(kind=c_int) :: fourthVertexID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_quad + + subroutine precicef_set_mesh_quads(meshName, size, ids, meshNameLength) & + & bind(c, name='precicef_set_mesh_quads_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_quads + + subroutine precicef_set_tetrahedron(meshName, firstVertexID, secondVertexID, & + & thirdVertexID, fourthVertexID, & + & meshNameLength) & + & bind(c, name='precicef_set_tetrahedron_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: firstVertexID + integer(kind=c_int) :: secondVertexID + integer(kind=c_int) :: thirdVertexID + integer(kind=c_int) :: fourthVertexID + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_tetrahedron + + subroutine precicef_set_mesh_tetrahedra(meshName, size, ids, meshNameLength) & + & bind(c, name='precicef_set_mesh_tetrahedra_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_tetrahedra + + subroutine precicef_requires_initial_data(isRequired) & + & bind(c, name='precicef_requires_initial_data_') + + use, intrinsic :: iso_c_binding + integer(kind=c_int) :: isRequired + end subroutine precicef_requires_initial_data + + subroutine precicef_write_data(meshName, dataName, size, ids, & + & values, meshNameLength, dataNameLength) & + & bind(c, name='precicef_write_data_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: values(*) + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_write_data + + subroutine precicef_read_data(meshName, dataName, size, ids, & + & relativeReadTime, values, meshNameLength, & + & dataNameLength) & + & bind(c, name='precicef_read_data_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: relativeReadTime + real(kind=c_double) :: values(*) + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_read_data + + subroutine precicef_set_mesh_access_region(meshName, boundingBox, & + & meshNameLength) & + & bind(c, name='precicef_set_mesh_access_region_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + real(kind=c_double) :: boundingBox(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_set_mesh_access_region + + subroutine precicef_get_mesh_vertex_ids_and_coordinates(meshName, size, & + & ids, coordinates, & + & meshNameLength) & + & bind(c, name='precicef_get_mesh_vertex_ids_and_coordinates_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: coordinates(*) + integer(kind=c_int), value :: meshNameLength + end subroutine precicef_get_mesh_vertex_ids_and_coordinates + + subroutine precicef_requires_gradient_data_for(meshName, dataName, & + & required, meshNameLength, & + & dataNameLength) & + & bind(c, name='precicef_requires_gradient_data_for_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: required + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_requires_gradient_data_for + + subroutine precicef_write_gradient_data(meshName, dataName, size, ids, & + & gradients, meshNameLength, & + & dataNameLength) & + & bind(c, name='precicef_write_gradient_data_') + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: meshName + character(kind=c_char), dimension(*) :: dataName + integer(kind=c_int) :: size + integer(kind=c_int) :: ids(*) + real(kind=c_double) :: gradients(*) + integer(kind=c_int), value :: meshNameLength + integer(kind=c_int), value :: dataNameLength + end subroutine precicef_write_gradient_data + + subroutine precicef_get_version_information(versionInfo, lengthVersionInfo) & + & bind(c, name="precicef_get_version_information_") + + use, intrinsic :: iso_c_binding + character(kind=c_char), dimension(*) :: versionInfo + integer(kind=c_int), value :: lengthVersionInfo + end subroutine precicef_get_version_information + + end interface + +end module precice From 05e57516b7bb869c9e5f67951d133de85f04d87f Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Tue, 7 Jan 2025 16:54:06 +0100 Subject: [PATCH 21/24] Add Fortran module and update precicef_ calls with length args --- elastic-tube-1d/fluid-fortran/CMakeLists.txt | 1 + .../fluid-fortran/src/FluidSolver.f90 | 18 ++++++++++------- elastic-tube-1d/solid-fortran/CMakeLists.txt | 1 + .../solid-fortran/src/SolidSolver.f90 | 20 +++++++++++++------ 4 files changed, 27 insertions(+), 13 deletions(-) diff --git a/elastic-tube-1d/fluid-fortran/CMakeLists.txt b/elastic-tube-1d/fluid-fortran/CMakeLists.txt index 1978e7d9a..bdb070435 100644 --- a/elastic-tube-1d/fluid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/fluid-fortran/CMakeLists.txt @@ -22,6 +22,7 @@ add_executable(FluidSolver src/FluidComputeSolution.f90 src/utilities.f90 src/FluidSolver.f90 + src/precice.f90 ) target_link_libraries(FluidSolver PRIVATE precice::precice) diff --git a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 index f59b30df4..8aa4833cb 100644 --- a/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 +++ b/elastic-tube-1d/fluid-fortran/src/FluidSolver.f90 @@ -1,4 +1,5 @@ program FluidSolver + use precice use FluidComputeSolution, only: fluid_compute_solution use Utilities, only: write_vtk implicit none @@ -39,7 +40,7 @@ program FluidSolver ! Configure precice rank = 0 commsize = 1 - call precicef_create(solverName, configFileName, rank, commsize) + call precicef_create(solverName, configFileName, rank, commsize, len_trim(solverName), len_trim(configFileName)) write(*, *) "preCICE configured..." ! Define mesh and data names @@ -53,7 +54,7 @@ program FluidSolver l = 10.0_dp ! Get mesh dimensions - call precicef_get_mesh_dimensions(meshName, dimensions) + call precicef_get_mesh_dimensions(meshName, dimensions, len_trim(meshName)) ! Allocate arrays allocate(vertexIDs(chunkLength)) @@ -100,7 +101,7 @@ program FluidSolver vertexIDs(i) = i - 1 end do - call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) + call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs, len_trim(meshName)) ! Check if Initial Data is Required and Write if Necessary call precicef_requires_initial_data(bool) @@ -112,7 +113,8 @@ program FluidSolver call precicef_initialize() ! read initial cross-Section length - call precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, 0.0d0, crossSectionLength) + call precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, 0.0d0, crossSectionLength, & + len_trim(meshName), len_trim(crossSectionLengthName)) ! Copy current cross-Section length to old array crossSectionLength_old = crossSectionLength @@ -147,13 +149,15 @@ program FluidSolver velocity, pressure, & ! resulting velocity pressure info) - call precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure) - + call precicef_write_data(meshName, pressureName, chunkLength, vertexIDs, pressure, & + len_trim(meshName), len_trim(pressureName)) + call precicef_advance(dt) call precicef_get_max_time_step_size(dt) - call precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, dt, crossSectionLength) + call precicef_read_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, dt, crossSectionLength, & + len_trim(meshName), len_trim(crossSectionLengthName)) call precicef_requires_reading_checkpoint(bool) if (bool .eq. 1) then diff --git a/elastic-tube-1d/solid-fortran/CMakeLists.txt b/elastic-tube-1d/solid-fortran/CMakeLists.txt index c33915351..7f57ccca4 100644 --- a/elastic-tube-1d/solid-fortran/CMakeLists.txt +++ b/elastic-tube-1d/solid-fortran/CMakeLists.txt @@ -20,6 +20,7 @@ find_package(precice 3 REQUIRED CONFIG) add_executable(SolidSolver src/SolidComputeSolution.f90 src/SolidSolver.f90 + src/precice.f90 ) target_link_libraries(SolidSolver PRIVATE precice::precice) diff --git a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 index 9de1b1d63..86be2925f 100644 --- a/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 +++ b/elastic-tube-1d/solid-fortran/src/SolidSolver.f90 @@ -1,4 +1,5 @@ program SolidSolver + use precice use SolidComputeSolution, only: solid_compute_solution implicit none @@ -39,10 +40,12 @@ program SolidSolver rank = 0 commsize = 1 - call precicef_create(solverName, configFileName, rank, commsize) + call precicef_create(solverName, configFileName, rank, commsize, & + len_trim(solverName), len_trim(configFileName)) write(*, *) 'preCICE configured...' - call precicef_get_mesh_dimensions(meshName, dimensions) + call precicef_get_mesh_dimensions(meshName, dimensions, & + len_trim(meshName)) ! Allocate arrays allocate(pressure(chunkLength)) @@ -53,18 +56,21 @@ program SolidSolver pressure = 0.0_dp crossSectionLength = 1.0_dp dx = tubeLength / real(domainSize, dp) + do i = 1, chunkLength grid((i - 1)*dimensions + 1) = dx * real(i - 1, dp) ! x-coordinate grid((i - 1)*dimensions + 2) = 0.0_dp ! y-coordinate vertexIDs(i) = i - 1 ! 0-based indexing here end do - call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs) + call precicef_set_vertices(meshName, chunkLength, grid, vertexIDs, & + len_trim(meshName)) ! Check if initial data is required and write if necessary call precicef_requires_initial_data(bool) if (bool == 1) then - call precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) + call precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, & + crossSectionLength, len_trim(meshName), len_trim(crossSectionLengthName)) end if write (*, *) 'Initialize preCICE...' @@ -81,11 +87,13 @@ program SolidSolver call precicef_get_max_time_step_size(dt) - call precicef_read_data(meshName, pressureName, chunkLength, vertexIDs, dt, pressure) + call precicef_read_data(meshName, pressureName, chunkLength, vertexIDs, dt, pressure, & + len_trim(meshName), len_trim(pressureName)) call solid_compute_solution(chunkLength, pressure, crossSectionLength) - call precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, crossSectionLength) + call precicef_write_data(meshName, crossSectionLengthName, chunkLength, vertexIDs, & + crossSectionLength, len_trim(meshName), len_trim(crossSectionLengthName)) call precicef_advance(dt) From 597f442223c4f7a8b28127af9e35a4e7c0bd36a6 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Wed, 8 Jan 2025 11:16:41 +0100 Subject: [PATCH 22/24] Add plotting for fluid fortran diameter --- elastic-tube-1d/plot-diameter.sh | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/elastic-tube-1d/plot-diameter.sh b/elastic-tube-1d/plot-diameter.sh index ae5f84d97..059eba585 100755 --- a/elastic-tube-1d/plot-diameter.sh +++ b/elastic-tube-1d/plot-diameter.sh @@ -23,3 +23,10 @@ if [ -n "$(ls -A ./fluid-rust/output/*.vtk 2>/dev/null)" ]; then else echo "No results to plot from fluid-python." fi + +# Plot diameter from fluid-fortran +if [ -n "$(ls -A ./fluid-fortran/output/*.vtk 2>/dev/null)" ]; then + python3 plot-vtk.py diameter fluid-fortran/output/out_fluid_ & +else + echo "No results to plot from fluid-fortran." +fi \ No newline at end of file From ab78cb6a84dc193cd8cd03c5d7a8079d4ffaf1f0 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Wed, 8 Jan 2025 11:54:28 +0100 Subject: [PATCH 23/24] Delete file (now fetched dynamically) --- elastic-tube-1d/fluid-fortran/src/precice.f90 | 315 ------------------ elastic-tube-1d/solid-fortran/src/precice.f90 | 315 ------------------ 2 files changed, 630 deletions(-) delete mode 100644 elastic-tube-1d/fluid-fortran/src/precice.f90 delete mode 100644 elastic-tube-1d/solid-fortran/src/precice.f90 diff --git a/elastic-tube-1d/fluid-fortran/src/precice.f90 b/elastic-tube-1d/fluid-fortran/src/precice.f90 deleted file mode 100644 index 13f880e2b..000000000 --- a/elastic-tube-1d/fluid-fortran/src/precice.f90 +++ /dev/null @@ -1,315 +0,0 @@ -module precice - use, intrinsic :: iso_c_binding - implicit none - - interface - - subroutine precicef_create(participantName, configFileName, & - & solverProcessIndex, solverProcessSize, & - & participantNameLength, configFileNameLength) & - & bind(c, name='precicef_create_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: participantName - character(kind=c_char), dimension(*) :: configFileName - integer(kind=c_int) :: solverProcessIndex - integer(kind=c_int) :: solverProcessSize - integer(kind=c_int), value :: participantNameLength - integer(kind=c_int), value :: configFileNameLength - end subroutine precicef_create - - subroutine precicef_initialize() & - & bind(c, name='precicef_initialize_') - - use, intrinsic :: iso_c_binding - end subroutine precicef_initialize - - subroutine precicef_advance(timestepLengthLimit) & - & bind(c, name='precicef_advance_') - - use, intrinsic :: iso_c_binding - real(kind=c_double) :: timestepLengthLimit - end subroutine precicef_advance - - subroutine precicef_finalize() & - & bind(c, name='precicef_finalize_') - - use, intrinsic :: iso_c_binding - end subroutine precicef_finalize - - subroutine precicef_requires_reading_checkpoint(isRequired) & - & bind(c, name='precicef_requires_reading_checkpoint_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isRequired - end subroutine precicef_requires_reading_checkpoint - - subroutine precicef_requires_writing_checkpoint(isRequired) & - & bind(c, name='precicef_requires_writing_checkpoint_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isRequired - end subroutine precicef_requires_writing_checkpoint - - subroutine precicef_get_mesh_dimensions(meshName, dimensions, & - & meshNameLength) & - & bind(c, name='precicef_get_mesh_dimensions_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: dimensions - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_get_mesh_dimensions - - subroutine precicef_get_data_dimensions(meshName, dataName, & - & dimensions, meshNameLength, dataNameLength) & - & bind(c, name='precicef_get_data_dimensions_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: dimensions - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_get_data_dimensions - - subroutine precicef_is_coupling_ongoing(isOngoing) & - & bind(c, name='precicef_is_coupling_ongoing_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isOngoing - end subroutine precicef_is_coupling_ongoing - - subroutine precicef_is_time_window_complete(isComplete) & - & bind(c, name='precicef_is_time_window_complete_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isComplete - end subroutine precicef_is_time_window_complete - - subroutine precicef_get_max_time_step_size(maxTimeStepSize) & - & bind(c, name='precicef_get_max_time_step_size_') - - use, intrinsic :: iso_c_binding - real(kind=c_double) :: maxTimeStepSize - end subroutine precicef_get_max_time_step_size - - subroutine precicef_requires_mesh_connectivity_for(meshName, required, meshNameLength) & - & bind(c, name='precicef_requires_mesh_connectivity_for_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: required - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_requires_mesh_connectivity_for - - subroutine precicef_set_vertex(meshName, position, vertexID, meshNameLength) & - & bind(c, name='precicef_set_vertex_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - real(kind=c_double) :: coordinates(3) - integer(kind=c_int) :: id - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_vertex - - subroutine precicef_get_mesh_vertex_size(meshName, meshSize, meshNameLength) & - & bind(c, name='precicef_get_mesh_vertex_size_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: meshSize - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_get_mesh_vertex_size - - subroutine precicef_set_vertices(meshName, size, coordinates, ids, meshNameLength) & - & bind(c, name='precicef_set_vertices_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - real(kind=c_double) :: coordinates(*) - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_vertices - - subroutine precicef_set_edge(meshName, firstVertexID, secondVertexID, & - & meshNameLength) & - & bind(c, name='precicef_set_edge_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstVertexID - integer(kind=c_int) :: secondVertexID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_edge - - subroutine precicef_set_mesh_edges(meshName, size, ids, meshNameLength) & - & bind(c, name='precicef_set_mesh_edges_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_edges - - subroutine precicef_set_triangle(meshName, firstEdgeID, secondEdgeID, & - & thirdEdgeID, meshNameLength) & - & bind(c, name='precicef_set_triangle_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstEdgeID - integer(kind=c_int) :: secondEdgeID - integer(kind=c_int) :: thirdEdgeID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_triangle - - subroutine precicef_set_quad(meshName, firstVertexID, secondVertexID, & - & thirdVertexID, fourthVertexID, & - & meshNameLength ) & - & bind(c, name='precicef_set_quad_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstVertexID - integer(kind=c_int) :: secondVertexID - integer(kind=c_int) :: thirdVertexID - integer(kind=c_int) :: fourthVertexID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_quad - - subroutine precicef_set_mesh_quads(meshName, size, ids, meshNameLength) & - & bind(c, name='precicef_set_mesh_quads_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_quads - - subroutine precicef_set_tetrahedron(meshName, firstVertexID, secondVertexID, & - & thirdVertexID, fourthVertexID, & - & meshNameLength) & - & bind(c, name='precicef_set_tetrahedron_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstVertexID - integer(kind=c_int) :: secondVertexID - integer(kind=c_int) :: thirdVertexID - integer(kind=c_int) :: fourthVertexID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_tetrahedron - - subroutine precicef_set_mesh_tetrahedra(meshName, size, ids, meshNameLength) & - & bind(c, name='precicef_set_mesh_tetrahedra_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_tetrahedra - - subroutine precicef_requires_initial_data(isRequired) & - & bind(c, name='precicef_requires_initial_data_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isRequired - end subroutine precicef_requires_initial_data - - subroutine precicef_write_data(meshName, dataName, size, ids, & - & values, meshNameLength, dataNameLength) & - & bind(c, name='precicef_write_data_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: values(*) - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_write_data - - subroutine precicef_read_data(meshName, dataName, size, ids, & - & relativeReadTime, values, meshNameLength, & - & dataNameLength) & - & bind(c, name='precicef_read_data_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: relativeReadTime - real(kind=c_double) :: values(*) - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_read_data - - subroutine precicef_set_mesh_access_region(meshName, boundingBox, & - & meshNameLength) & - & bind(c, name='precicef_set_mesh_access_region_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - real(kind=c_double) :: boundingBox(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_access_region - - subroutine precicef_get_mesh_vertex_ids_and_coordinates(meshName, size, & - & ids, coordinates, & - & meshNameLength) & - & bind(c, name='precicef_get_mesh_vertex_ids_and_coordinates_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: coordinates(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_get_mesh_vertex_ids_and_coordinates - - subroutine precicef_requires_gradient_data_for(meshName, dataName, & - & required, meshNameLength, & - & dataNameLength) & - & bind(c, name='precicef_requires_gradient_data_for_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: required - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_requires_gradient_data_for - - subroutine precicef_write_gradient_data(meshName, dataName, size, ids, & - & gradients, meshNameLength, & - & dataNameLength) & - & bind(c, name='precicef_write_gradient_data_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: gradients(*) - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_write_gradient_data - - subroutine precicef_get_version_information(versionInfo, lengthVersionInfo) & - & bind(c, name="precicef_get_version_information_") - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: versionInfo - integer(kind=c_int), value :: lengthVersionInfo - end subroutine precicef_get_version_information - - end interface - -end module precice diff --git a/elastic-tube-1d/solid-fortran/src/precice.f90 b/elastic-tube-1d/solid-fortran/src/precice.f90 deleted file mode 100644 index 13f880e2b..000000000 --- a/elastic-tube-1d/solid-fortran/src/precice.f90 +++ /dev/null @@ -1,315 +0,0 @@ -module precice - use, intrinsic :: iso_c_binding - implicit none - - interface - - subroutine precicef_create(participantName, configFileName, & - & solverProcessIndex, solverProcessSize, & - & participantNameLength, configFileNameLength) & - & bind(c, name='precicef_create_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: participantName - character(kind=c_char), dimension(*) :: configFileName - integer(kind=c_int) :: solverProcessIndex - integer(kind=c_int) :: solverProcessSize - integer(kind=c_int), value :: participantNameLength - integer(kind=c_int), value :: configFileNameLength - end subroutine precicef_create - - subroutine precicef_initialize() & - & bind(c, name='precicef_initialize_') - - use, intrinsic :: iso_c_binding - end subroutine precicef_initialize - - subroutine precicef_advance(timestepLengthLimit) & - & bind(c, name='precicef_advance_') - - use, intrinsic :: iso_c_binding - real(kind=c_double) :: timestepLengthLimit - end subroutine precicef_advance - - subroutine precicef_finalize() & - & bind(c, name='precicef_finalize_') - - use, intrinsic :: iso_c_binding - end subroutine precicef_finalize - - subroutine precicef_requires_reading_checkpoint(isRequired) & - & bind(c, name='precicef_requires_reading_checkpoint_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isRequired - end subroutine precicef_requires_reading_checkpoint - - subroutine precicef_requires_writing_checkpoint(isRequired) & - & bind(c, name='precicef_requires_writing_checkpoint_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isRequired - end subroutine precicef_requires_writing_checkpoint - - subroutine precicef_get_mesh_dimensions(meshName, dimensions, & - & meshNameLength) & - & bind(c, name='precicef_get_mesh_dimensions_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: dimensions - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_get_mesh_dimensions - - subroutine precicef_get_data_dimensions(meshName, dataName, & - & dimensions, meshNameLength, dataNameLength) & - & bind(c, name='precicef_get_data_dimensions_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: dimensions - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_get_data_dimensions - - subroutine precicef_is_coupling_ongoing(isOngoing) & - & bind(c, name='precicef_is_coupling_ongoing_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isOngoing - end subroutine precicef_is_coupling_ongoing - - subroutine precicef_is_time_window_complete(isComplete) & - & bind(c, name='precicef_is_time_window_complete_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isComplete - end subroutine precicef_is_time_window_complete - - subroutine precicef_get_max_time_step_size(maxTimeStepSize) & - & bind(c, name='precicef_get_max_time_step_size_') - - use, intrinsic :: iso_c_binding - real(kind=c_double) :: maxTimeStepSize - end subroutine precicef_get_max_time_step_size - - subroutine precicef_requires_mesh_connectivity_for(meshName, required, meshNameLength) & - & bind(c, name='precicef_requires_mesh_connectivity_for_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: required - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_requires_mesh_connectivity_for - - subroutine precicef_set_vertex(meshName, position, vertexID, meshNameLength) & - & bind(c, name='precicef_set_vertex_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - real(kind=c_double) :: coordinates(3) - integer(kind=c_int) :: id - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_vertex - - subroutine precicef_get_mesh_vertex_size(meshName, meshSize, meshNameLength) & - & bind(c, name='precicef_get_mesh_vertex_size_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: meshSize - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_get_mesh_vertex_size - - subroutine precicef_set_vertices(meshName, size, coordinates, ids, meshNameLength) & - & bind(c, name='precicef_set_vertices_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - real(kind=c_double) :: coordinates(*) - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_vertices - - subroutine precicef_set_edge(meshName, firstVertexID, secondVertexID, & - & meshNameLength) & - & bind(c, name='precicef_set_edge_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstVertexID - integer(kind=c_int) :: secondVertexID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_edge - - subroutine precicef_set_mesh_edges(meshName, size, ids, meshNameLength) & - & bind(c, name='precicef_set_mesh_edges_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_edges - - subroutine precicef_set_triangle(meshName, firstEdgeID, secondEdgeID, & - & thirdEdgeID, meshNameLength) & - & bind(c, name='precicef_set_triangle_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstEdgeID - integer(kind=c_int) :: secondEdgeID - integer(kind=c_int) :: thirdEdgeID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_triangle - - subroutine precicef_set_quad(meshName, firstVertexID, secondVertexID, & - & thirdVertexID, fourthVertexID, & - & meshNameLength ) & - & bind(c, name='precicef_set_quad_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstVertexID - integer(kind=c_int) :: secondVertexID - integer(kind=c_int) :: thirdVertexID - integer(kind=c_int) :: fourthVertexID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_quad - - subroutine precicef_set_mesh_quads(meshName, size, ids, meshNameLength) & - & bind(c, name='precicef_set_mesh_quads_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_quads - - subroutine precicef_set_tetrahedron(meshName, firstVertexID, secondVertexID, & - & thirdVertexID, fourthVertexID, & - & meshNameLength) & - & bind(c, name='precicef_set_tetrahedron_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: firstVertexID - integer(kind=c_int) :: secondVertexID - integer(kind=c_int) :: thirdVertexID - integer(kind=c_int) :: fourthVertexID - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_tetrahedron - - subroutine precicef_set_mesh_tetrahedra(meshName, size, ids, meshNameLength) & - & bind(c, name='precicef_set_mesh_tetrahedra_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_tetrahedra - - subroutine precicef_requires_initial_data(isRequired) & - & bind(c, name='precicef_requires_initial_data_') - - use, intrinsic :: iso_c_binding - integer(kind=c_int) :: isRequired - end subroutine precicef_requires_initial_data - - subroutine precicef_write_data(meshName, dataName, size, ids, & - & values, meshNameLength, dataNameLength) & - & bind(c, name='precicef_write_data_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: values(*) - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_write_data - - subroutine precicef_read_data(meshName, dataName, size, ids, & - & relativeReadTime, values, meshNameLength, & - & dataNameLength) & - & bind(c, name='precicef_read_data_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: relativeReadTime - real(kind=c_double) :: values(*) - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_read_data - - subroutine precicef_set_mesh_access_region(meshName, boundingBox, & - & meshNameLength) & - & bind(c, name='precicef_set_mesh_access_region_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - real(kind=c_double) :: boundingBox(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_set_mesh_access_region - - subroutine precicef_get_mesh_vertex_ids_and_coordinates(meshName, size, & - & ids, coordinates, & - & meshNameLength) & - & bind(c, name='precicef_get_mesh_vertex_ids_and_coordinates_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: coordinates(*) - integer(kind=c_int), value :: meshNameLength - end subroutine precicef_get_mesh_vertex_ids_and_coordinates - - subroutine precicef_requires_gradient_data_for(meshName, dataName, & - & required, meshNameLength, & - & dataNameLength) & - & bind(c, name='precicef_requires_gradient_data_for_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: required - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_requires_gradient_data_for - - subroutine precicef_write_gradient_data(meshName, dataName, size, ids, & - & gradients, meshNameLength, & - & dataNameLength) & - & bind(c, name='precicef_write_gradient_data_') - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: meshName - character(kind=c_char), dimension(*) :: dataName - integer(kind=c_int) :: size - integer(kind=c_int) :: ids(*) - real(kind=c_double) :: gradients(*) - integer(kind=c_int), value :: meshNameLength - integer(kind=c_int), value :: dataNameLength - end subroutine precicef_write_gradient_data - - subroutine precicef_get_version_information(versionInfo, lengthVersionInfo) & - & bind(c, name="precicef_get_version_information_") - - use, intrinsic :: iso_c_binding - character(kind=c_char), dimension(*) :: versionInfo - integer(kind=c_int), value :: lengthVersionInfo - end subroutine precicef_get_version_information - - end interface - -end module precice From 1b754050dcf1ab9481d71e084a1ebcbdc59181a4 Mon Sep 17 00:00:00 2001 From: YonatanGM Date: Wed, 8 Jan 2025 12:00:05 +0100 Subject: [PATCH 24/24] Get precice.f90 from the fortran module repo in run script --- elastic-tube-1d/fluid-fortran/run.sh | 5 +++++ elastic-tube-1d/solid-fortran/run.sh | 5 +++++ 2 files changed, 10 insertions(+) diff --git a/elastic-tube-1d/fluid-fortran/run.sh b/elastic-tube-1d/fluid-fortran/run.sh index 2a0de7786..1f307ae6f 100755 --- a/elastic-tube-1d/fluid-fortran/run.sh +++ b/elastic-tube-1d/fluid-fortran/run.sh @@ -4,6 +4,11 @@ set -e -u . ../../tools/log.sh exec > >(tee --append "$LOGFILE") 2>&1 +if [ ! -f src/precice.f90 ]; then + echo "Fetching precice.f90 (Module for Fortran bindings of preCICE)..." + curl -o src/precice.f90 https://raw.githubusercontent.com/precice/fortran-module/master/precice.f90 +fi + if [ ! -d build ]; then mkdir build cmake -S . -B build diff --git a/elastic-tube-1d/solid-fortran/run.sh b/elastic-tube-1d/solid-fortran/run.sh index 5263549b9..337bd8c1f 100755 --- a/elastic-tube-1d/solid-fortran/run.sh +++ b/elastic-tube-1d/solid-fortran/run.sh @@ -4,6 +4,11 @@ set -e -u . ../../tools/log.sh || true exec > >(tee --append "$LOGFILE") 2>&1 || true +if [ ! -f src/precice.f90 ]; then + echo "Fetching precice.f90 (Module for Fortran bindings of preCICE)..." + curl -o src/precice.f90 https://raw.githubusercontent.com/precice/fortran-module/master/precice.f90 +fi + if [ ! -d build ]; then mkdir build cd build