Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 9 additions & 2 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
CMP = gfortran# ifort,ifx,gfortran

BUILD ?=
PRECISION ?= single

#######CMP settings###########
ifeq ($(CMP),ifort)
Expand All @@ -14,9 +15,9 @@ FFLAGS = -diag-disable=10448 -fpp -O3 -free -qopenmp -heap-arrays #-real-size 32
else ifeq ($(CMP),ifx)
FC = ifx
FFLAGS = -diag-disable=10448 -fpp -O3 -free -qopenmp -heap-arrays #-real-size 32 -double-size 64
else ifeq ($(CMP),gfortran)
else ifeq ($(filter $(CMP),gfortran gcc),$(CMP))
FC = gfortran
FFLAGS = -O3
FFLAGS = -cpp -O3
ifeq ($(BUILD),debug)
FFLAGS = -cpp -g -O0
FFLAGS += -ffpe-trap=invalid,zero -fbacktrace -Wall -Wextra -pedantic -Warray-bounds -fbacktrace -fbounds-check
Expand All @@ -26,6 +27,12 @@ endif
CC = cc
CPP = c++

ifeq ($(PRECISION),single)
FFLAGS += -DFSILBM_SINGLE_PRECISION
else ifneq ($(PRECISION),double)
$(error PRECISION must be either single or double)
endif

SRCDIR = ./src

### List of files for the main code
Expand Down
34 changes: 20 additions & 14 deletions src/ConstParams.f90
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,12 @@
! Copyright (C) 2025-2026 Ankang Gao and contributors

module ConstParams
integer, parameter :: sp = selected_real_kind(6, 37)
#ifdef FSILBM_SINGLE_PRECISION
integer, parameter :: rp = sp
#else
integer, parameter :: rp = selected_real_kind(15, 307)
#endif
! LBM module
! D3Q19model
integer, parameter:: SpaceDim = 3, lbmDim = 18
Expand All @@ -19,28 +25,28 @@ module ConstParams
integer, parameter:: positivedirs(1:lbmDim/2) = [1, 3, 5, 7, 8, 11, 12, 15, 16]
integer, parameter:: negativedirs(1:lbmDim/2) = [2, 4, 6,10, 9, 14, 13, 18, 17]
! Weights
real(8), parameter:: wt(0:lbmDim) = [&
1.0d0/3.0d0,1.0d0/18.0d0,1.0d0/18.0d0,1.0d0/18.0d0,1.0d0/18.0d0,1.0d0/18.0d0,1.0d0/18.0d0, &
1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0, &
1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0 ]
real(rp), parameter:: wt(0:lbmDim) = [&
1.0e0_rp/3.0e0_rp,1.0e0_rp/18.0e0_rp,1.0e0_rp/18.0e0_rp,1.0e0_rp/18.0e0_rp,1.0e0_rp/18.0e0_rp,1.0e0_rp/18.0e0_rp,1.0e0_rp/18.0e0_rp, &
1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp, &
1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp,1.0e0_rp/36.0e0_rp ]
! MRT model Matrix
!real(8),parameter::s0=1.0d0,s1=1.0d0,s2=1.0d0,s4=1.0d0,s10=1.0d0,s16=1.0d0
real(8), parameter::s0=0.0d0,s1=1.19d0,s2=1.4d0,s4=1.2d0,s10=1.4d0,s16=1.98d0
!real(rp),parameter::s0=1.0e0_rp,s1=1.0e0_rp,s2=1.0e0_rp,s4=1.0e0_rp,s10=1.0e0_rp,s16=1.0e0_rp
real(rp), parameter::s0=0.0e0_rp,s1=1.19e0_rp,s2=1.4e0_rp,s4=1.2e0_rp,s10=1.4e0_rp,s16=1.98e0_rp

! MODULE BoundCondParams
integer, parameter:: BCEq_DirecletU = 101,BCnEq_DirecletU = 102,BCorder1_Extrapolate = 103,BCorder2_Extrapolate = 104
!given balance function,unbalanced extrapolation,1st order extrapolate,2nd order extrapolate
integer, parameter:: BCstationary_Wall = 201, BCmoving_Wall = 202, BCstationary_Wall_halfway = 203, BCmoving_Wall_halfway = 204
integer, parameter:: BCPeriodic = 301,BCSymmetric = 302,BCfluid = 0,BCfluid_father = 1

real(8), parameter:: Pi = 3.141592653589793d0,eps = 1.0d-5,MachineTolerace = 1.0d-12
real(rp), parameter:: Pi = 3.141592653589793e0_rp,eps = 1.0e-5_rp,MachineTolerace = 1.0e-12_rp
integer, parameter:: DOFDim = 6

real(8), parameter:: Cs2 = 1.d0/3.0d0
real(8), parameter:: Csmag = 0.17d0
! real(8), parameter:: CsmagConst = 16.d0 * dsqrt(2.d0) / (3.d0 * Pi * Pi)
real(8), parameter:: CsmagConst = 2.0d0 * Csmag * Csmag * dsqrt(2.d0) * 9.d0
real(8), parameter:: CWALE = 0.50d0
real(8), parameter:: CWALEConst = CWALE*CWALE
real(8), parameter:: CvremConst = 2.5d0*Csmag*Csmag
real(rp), parameter:: Cs2 = 1.e0_rp/3.0e0_rp
real(rp), parameter:: Csmag = 0.17e0_rp
! real(rp), parameter:: CsmagConst = 16.e0_rp * sqrt(2.e0_rp) / (3.e0_rp * Pi * Pi)
real(rp), parameter:: CsmagConst = 2.0e0_rp * Csmag * Csmag * sqrt(2.e0_rp) * 9.e0_rp
real(rp), parameter:: CWALE = 0.50e0_rp
real(rp), parameter:: CWALEConst = CWALE*CWALE
real(rp), parameter:: CvremConst = 2.5e0_rp*Csmag*Csmag
end module ConstParams
25 changes: 13 additions & 12 deletions src/FlowCondition.f90
Original file line number Diff line number Diff line change
Expand Up @@ -4,25 +4,26 @@
! Copyright (C) 2025-2026 Ankang Gao and contributors

module FlowCondition
use ConstParams, only: rp
implicit none
private
public :: FlowCondType,flow
public :: read_flow_conditions,read_probe_params,write_information_titles,write_fluid_information
type :: FlowCondType
integer :: isConCmpt,numsubstep,npsize
real(8) :: timeSimTotal,timeContiDelta,timeWriteBegin,timeWriteEnd,timeFlowDelta,timeBodyDelta,timeInfoDelta
real(8) :: Re,denIn,nu,Mu,dtolLBM
real(rp) :: timeSimTotal,timeContiDelta,timeWriteBegin,timeWriteEnd,timeFlowDelta,timeBodyDelta,timeInfoDelta
real(rp) :: Re,denIn,nu,Mu,dtolLBM
integer :: LrefType,TrefType,UrefType,ntolLBM
integer :: velocityKind,interpolateScheme
real(8) :: uvwIn(1:3),shearRateIn(1:3)
real(8) :: volumeForceIn(1:3),volumeForceAmp,volumeForceFreq,volumeForcePhi
real(8) :: Uref,Lref,Tref
real(8) :: Aref,Fref,Eref,Pref
real(8) :: Asfac,Lchod,Lspan,AR
real(rp) :: uvwIn(1:3),shearRateIn(1:3)
real(rp) :: volumeForceIn(1:3),volumeForceAmp,volumeForceFreq,volumeForcePhi
real(rp) :: Uref,Lref,Tref
real(rp) :: Aref,Fref,Eref,Pref
real(rp) :: Asfac,Lchod,Lspan,AR
integer :: fluidProbingNum,inWhichBlock,solidProbingNum
integer, allocatable :: solidProbingNode(:)
real(8), allocatable :: fluidProbingCoords(:,:)
real(8) :: AmplInitDist(1:3),waveInitDist
real(rp), allocatable :: fluidProbingCoords(:,:)
real(rp) :: AmplInitDist(1:3),waveInitDist
end type FlowCondType
type(FlowCondType) :: flow

Expand Down Expand Up @@ -71,7 +72,7 @@ SUBROUTINE read_flow_conditions(filename)
read(buffer,*) flow%interpolateScheme
close(111)
! flow%denIn is not 1
if(abs(flow%denIn-1.d0).gt.1e-6) then
if(abs(flow%denIn-1.e0_rp).gt.1e-6_rp) then
write(*,*) 'Warning, denIn is not 1, ', flow%denIn
endif
END SUBROUTINE
Expand Down Expand Up @@ -195,8 +196,8 @@ SUBROUTINE write_information_titles(nGroup)
SUBROUTINE write_fluid_information(time,dh,xmin,ymin,zmin,xDim,yDim,zDim,velocityIn)
implicit none
integer:: i,j,xDim,yDim,zDim
real(8):: time,dh,xmin,ymin,zmin,xmax,ymax,zmax
real(8):: velocityIn(zDim,yDim,xDim,1:3),velocityOut(1:3)
real(rp):: time,dh,xmin,ymin,zmin,xmax,ymax,zmax
real(rp):: velocityIn(zDim,yDim,xDim,1:3),velocityOut(1:3)
integer,parameter::probeLen=4
character (LEN=probeLen):: probeNum
! write fluid probing information
Expand Down
Loading