diff --git a/changelog-entries/607.md b/changelog-entries/607.md new file mode 100644 index 000000000..5f234dc48 --- /dev/null +++ b/changelog-entries/607.md @@ -0,0 +1 @@ +- Added Fortran solvers for the elastic-tube-1d tutorial using the [Fortran module](https://github.com/precice/fortran-module). \ No newline at end of file diff --git a/elastic-tube-1d/fluid-fortran-module/CMakeLists.txt b/elastic-tube-1d/fluid-fortran-module/CMakeLists.txt new file mode 100644 index 000000000..706fd822f --- /dev/null +++ b/elastic-tube-1d/fluid-fortran-module/CMakeLists.txt @@ -0,0 +1,28 @@ +cmake_minimum_required(VERSION 3.10) +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) + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS + "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +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) + +add_executable(FluidSolver + src/FluidComputeSolution.f90 + src/utilities.f90 + src/FluidSolver.f90 + src/precice.f90 +) + +target_link_libraries(FluidSolver PRIVATE precice::precice) +target_link_libraries(FluidSolver PRIVATE ${LAPACK_LIBRARIES} ${LAPACK_LINKER_FLAGS}) + diff --git a/elastic-tube-1d/fluid-fortran-module/clean.sh b/elastic-tube-1d/fluid-fortran-module/clean.sh new file mode 100755 index 000000000..8d5d37784 --- /dev/null +++ b/elastic-tube-1d/fluid-fortran-module/clean.sh @@ -0,0 +1,8 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +rm -rvf ./output/*.vtk +clean_precice_logs . +clean_case_logs . diff --git a/elastic-tube-1d/fluid-fortran-module/run.sh b/elastic-tube-1d/fluid-fortran-module/run.sh new file mode 100755 index 000000000..0b6d315c7 --- /dev/null +++ b/elastic-tube-1d/fluid-fortran-module/run.sh @@ -0,0 +1,22 @@ +#!/usr/bin/env bash +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 + cmake --build build +fi + +mkdir -p output + +./build/FluidSolver ../precice-config.xml + +close_log diff --git a/elastic-tube-1d/fluid-fortran-module/src/FluidComputeSolution.f90 b/elastic-tube-1d/fluid-fortran-module/src/FluidComputeSolution.f90 new file mode 100644 index 000000000..650fc79f6 --- /dev/null +++ b/elastic-tube-1d/fluid-fortran-module/src/FluidComputeSolution.f90 @@ -0,0 +1,182 @@ +module FluidComputeSolution + implicit none + integer, parameter :: dp = kind(1.0d0) + +contains + + subroutine fluid_compute_solution(velocity_old, pressure_old, & + crossSectionLength_old, crossSectionLength, t, N, kappa, tau, & + velocity, pressure, info) + + 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 + 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(dp) :: alpha, L, dx, velocity_in, tmp2, norm_1, norm_2, norm + + ! LAPACK Variables + integer :: nlhs, nrhs + real(dp), allocatable :: Res(:) + real(dp), 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); + + ! Output status from dgesv (0 = success, < 0 = invalid argument, > 0 = singular matrix) + info = 0 + + ! 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) + + call dgesv(nlhs, nrhs, LHS, nlhs, ipiv, Res, nlhs, info) + if (info /= 0) then + write(*, *) "Linear Solver not converged!, Info: ", info + 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 fluid_compute_solution +end module FluidComputeSolution diff --git a/elastic-tube-1d/fluid-fortran-module/src/FluidSolver.f90 b/elastic-tube-1d/fluid-fortran-module/src/FluidSolver.f90 new file mode 100644 index 000000000..f2b4485f5 --- /dev/null +++ b/elastic-tube-1d/fluid-fortran-module/src/FluidSolver.f90 @@ -0,0 +1,194 @@ +program FluidSolver + use precice + use FluidComputeSolution, only: fluid_compute_solution + use Utilities, only: write_vtk + implicit none + integer, parameter :: dp = kind(1.0d0) ! Double precision + + character(LEN=50) :: configFileName + character(LEN=50) :: solverName + character(LEN=50) :: meshName, pressureName, crossSectionLengthName + character(LEN=50) :: 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' + + ! Configure precice + rank = 0 + commsize = 1 + call precicef_create(solverName, configFileName, rank, commsize, len_trim(solverName), len_trim(configFileName)) + write(*, *) "preCICE configured..." + + ! Define mesh and data names + meshName = "Fluid-Nodes-Mesh" + pressureName = "Pressure" + crossSectionLengthName = "CrossSectionLength" + + domainSize = 100 + chunkLength = domainSize + 1 + kappa = 100.0_dp + l = 10.0_dp + + ! Get mesh dimensions + call precicef_get_mesh_dimensions(meshName, dimensions, len_trim(meshName)) + + ! 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) ! radius of the tube + a0 = r0**2 * PI ! cross-sectional area + u0 = 10.0_dp ! mean velocity + ampl = 3.0_dp ! amplitude of varying velocity + frequency = 10.0_dp ! frequency of variation + t_shift = 0.0_dp ! temporal shift of variation + p0 = 0.0_dp ! pressure at outlet + 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, 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 + + 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 + write (*, *) 'Fluid: Writing initial data' + end if + + write (*, *) "Initialize preCICE..." + call precicef_initialize() + + ! read initial cross-Section length + 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 + + ! 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 + + t = 0.0d0 + out_counter = 0 + + ! Main coupling loop + call precicef_is_coupling_ongoing(ongoing) + do while (ongoing /= 0) + ! checkpointing is required in implicit coupling + call precicef_requires_writing_checkpoint(bool) + if (bool .eq. 1) then + ! nothing + 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, & + dt, & ! tau + velocity, pressure, & ! resulting velocity pressure + info) + + 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, & + len_trim(meshName), len_trim(crossSectionLengthName)) + + call precicef_requires_reading_checkpoint(bool) + if (bool .eq. 1) then + ! not yet converged + else + ! converged, advance in time + 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 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-module/src/utilities.f90 b/elastic-tube-1d/fluid-fortran-module/src/utilities.f90 new file mode 100644 index 000000000..e844e0a64 --- /dev/null +++ b/elastic-tube-1d/fluid-fortran-module/src/utilities.f90 @@ -0,0 +1,74 @@ +module Utilities + implicit none + integer, parameter :: dp = kind(1.0D0) +contains + + subroutine write_vtk(t, iteration, filenamePrefix, nSlices, & + grid, velocity, pressure, diameter) + implicit none + + 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 + + write (filename, '(A,"_",I0,".vtk")') trim(filenamePrefix), iteration + 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 + print *, 'Error: Unable to open file ', trim(filename) + return + end if + + ! Write vtk headers + 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 + 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 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 + write (ioUnit, '(ES24.16)') pressure(i) + 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 + write (ioUnit, '(ES24.16)') diameter(i) + end do + write (ioUnit, '(A)') '' + + close (ioUnit) + + end subroutine write_vtk + +end module Utilities diff --git a/elastic-tube-1d/images/tutorials-elastic-tube-1d-all.png b/elastic-tube-1d/images/tutorials-elastic-tube-1d-all.png index 166b27d4e..0aaf61182 100644 Binary files a/elastic-tube-1d/images/tutorials-elastic-tube-1d-all.png and b/elastic-tube-1d/images/tutorials-elastic-tube-1d-all.png differ diff --git a/elastic-tube-1d/plot-all.sh b/elastic-tube-1d/plot-all.sh index 8ba00ee03..bb1396070 100755 --- a/elastic-tube-1d/plot-all.sh +++ b/elastic-tube-1d/plot-all.sh @@ -12,15 +12,18 @@ 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", \ + "fluid-fortran-module/precice-Fluid-watchpoint-Middle.log" using 1:4 with linespoints title "fortran-module" 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", \ + "fluid-fortran-module/precice-Fluid-watchpoint-Middle.log" using 1:5 with linespoints title "fortran-module" + + set xlabel 'time [s]' EOF diff --git a/elastic-tube-1d/plot-diameter.sh b/elastic-tube-1d/plot-diameter.sh index ae5f84d97..4f7c30d4e 100755 --- a/elastic-tube-1d/plot-diameter.sh +++ b/elastic-tube-1d/plot-diameter.sh @@ -23,3 +23,17 @@ 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 + +# Plot diameter from fluid-fortran-module +if [ -n "$(ls -A ./fluid-fortran-module/output/*.vtk 2>/dev/null)" ]; then + python3 plot-vtk.py diameter fluid-fortran-module/output/out_fluid_ & +else + echo "No results to plot from fluid-fortran-module." +fi \ No newline at end of file diff --git a/elastic-tube-1d/solid-fortran-module/CMakeLists.txt b/elastic-tube-1d/solid-fortran-module/CMakeLists.txt new file mode 100644 index 000000000..a1be8258f --- /dev/null +++ b/elastic-tube-1d/solid-fortran-module/CMakeLists.txt @@ -0,0 +1,24 @@ +cmake_minimum_required(VERSION 3.10) + +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) + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS + "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +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 + src/SolidComputeSolution.f90 + src/SolidSolver.f90 + src/precice.f90 +) + +target_link_libraries(SolidSolver PRIVATE precice::precice) diff --git a/elastic-tube-1d/solid-fortran-module/clean.sh b/elastic-tube-1d/solid-fortran-module/clean.sh new file mode 100755 index 000000000..ae7a54924 --- /dev/null +++ b/elastic-tube-1d/solid-fortran-module/clean.sh @@ -0,0 +1,7 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . diff --git a/elastic-tube-1d/solid-fortran-module/run.sh b/elastic-tube-1d/solid-fortran-module/run.sh new file mode 100755 index 000000000..337bd8c1f --- /dev/null +++ b/elastic-tube-1d/solid-fortran-module/run.sh @@ -0,0 +1,22 @@ +#!/usr/bin/env bash +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 + cmake .. + cmake --build . + cd .. +fi + +./build/SolidSolver ../precice-config.xml + +close_log diff --git a/elastic-tube-1d/solid-fortran-module/src/SolidComputeSolution.f90 b/elastic-tube-1d/solid-fortran-module/src/SolidComputeSolution.f90 new file mode 100644 index 000000000..28aeb72b4 --- /dev/null +++ b/elastic-tube-1d/solid-fortran-module/src/SolidComputeSolution.f90 @@ -0,0 +1,32 @@ +module SolidComputeSolution + implicit none + integer, parameter :: dp = kind(1.0d0) + +contains + + subroutine solid_compute_solution(chunkLength, pressure, crossSectionLength) + integer, intent(in) :: chunkLength + real(dp), intent(in) :: pressure(1:chunkLength) + real(dp), intent(inout) :: crossSectionLength(1:chunkLength) + + real(dp) :: pi, e, r0, c_mk, c_mk2 + real(dp) :: pressure0 + integer :: i + + ! 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.0_dp + + ! 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 solid_compute_solution + +end module SolidComputeSolution diff --git a/elastic-tube-1d/solid-fortran-module/src/SolidSolver.f90 b/elastic-tube-1d/solid-fortran-module/src/SolidSolver.f90 new file mode 100644 index 000000000..86be2925f --- /dev/null +++ b/elastic-tube-1d/solid-fortran-module/src/SolidSolver.f90 @@ -0,0 +1,117 @@ +program SolidSolver + use precice + use SolidComputeSolution, only: solid_compute_solution + implicit none + + integer, parameter :: dp = kind(1.0d0) + + 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 + 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.0_dp + + 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, & + len_trim(solverName), len_trim(configFileName)) + write(*, *) 'preCICE configured...' + + call precicef_get_mesh_dimensions(meshName, dimensions, & + len_trim(meshName)) + + ! 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 + + 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, len_trim(meshName), len_trim(crossSectionLengthName)) + 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 .eq. 1) then + ! Do nothing here + end if + + call precicef_get_max_time_step_size(dt) + + 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, len_trim(meshName), len_trim(crossSectionLengthName)) + + call precicef_advance(dt) + + call precicef_requires_reading_checkpoint(bool) + if (bool .eq. 1) then + ! nothing + end if + + call precicef_is_coupling_ongoing(ongoing) + end do + + write (*, *) 'Exiting SolidSolver' + + call precicef_finalize() + + deallocate(pressure) + deallocate(crossSectionLength) + deallocate(grid) + deallocate(vertexIDs) + +end program SolidSolver