Merge pull request #30 from chowland/ibm_rollback

IBM implementation + Direct interpolation of salinity to appropriate grid for buoyancy driving

Makefile

@@ -1,12 +1,12 @@
 # Choose the machine being used
-# Options: PC_GNU, PC_INTEL, (i)SNELLIUS, IRENE, MARENOSTRUM, SUPERMUC
+# Options: PC_GNU, PC_INTEL, (i)SNELLIUS, IRENE(_SKL/_ROME), MARENOSTRUM, SUPERMUC
 MACHINE=PC_GNU
 # Modules required for each HPC system as follows:
-# SNELLIUS: 2021 foss/2021a HDF5/1.10.7-gompi-2021a
-# iSNELLIUS: 2021 intel/2021a FFTW/3.3.9-intel-2021a HDF5/1.10.7-iimpi-2021a
+# SNELLIUS: 2022 foss/2022a HDF5/1.12.2-gompi-2022a
+# iSNELLIUS: 2022 intel/2022a FFTW/3.3.10-GCC-11.3.0 HDF5/1.12.2-iimpi-2021a
 # IRENE: flavor/hdf5/parallel hdf5 fftw3/gnu
 # MARENOSTRUM: fabric intel mkl impi hdf5 fftw szip
-# SUPERMUC: fftw hdf5 szip
+# SUPERMUC: fftw hdf5
 
 #=======================================================================
 #  Compiler options
@@ -29,22 +29,18 @@ ifeq ($(MACHINE),PC_INTEL)
 endif
 ifeq ($(MACHINE),iSNELLIUS)
 	FC = h5pfc -fpp -r8 -O3 -align array64byte -fma -ftz -fomit-frame-pointer
-	LDFLAGS = -lfftw3 -mkl=sequential
+	LDFLAGS = -lfftw3 -qmkl=sequential
 endif
 ifeq ($(MACHINE),SNELLIUS)
 	FC = h5pfc -cpp -fdefault-real-8 -fdefault-double-8 -w -fallow-argument-mismatch
-	FC += -O3 -march=znver1 -mtune=znver1 -mfma -mavx2 -m3dnow -fomit-frame-pointer
+	FC += -O2 -march=znver1 -mtune=znver1 -mfma -mavx2 -m3dnow -fomit-frame-pointer
 	BLAS_LIBS = -lscalapack -lopenblas -ldl
 	LDFLAGS = -lfftw3 $(BLAS_LIBS)
 endif
 ifeq ($(MACHINE),IRENE_SKL)
 	FC = h5pfc -fpp -r8 -O3 -mtune=skylake -xCORE-AVX512 -m64 -fPIC $(FFTW3_FFLAGS)
 	LDFLAGS = $(FFTW3_LDFLAGS) $(MKL_LDFLAGS) -ldl
 endif
-ifeq ($(MACHINE),IRENE_KNL)
-	FC = h5pfc -fpp -r8 -O3 -xMIC-AVX512 -fma -align array64byte -finline-functions $(FFTW3_FFLAGS)
-	LDFLAGS = $(FFTW3_LDFLAGS) $(MKL_LDFLAGS) -ldl
-endif
 ifeq ($(MACHINE),IRENE_ROME)
 	FC = h5pfc -fpp -r8 -O3 -mavx2 $(FFTW3_FFLAGS)
 	HDF5_LIBS = -lhdf5_fortran -lhdf5 -lz -ldl -lm
@@ -55,9 +51,8 @@ ifeq ($(MACHINE),MARENOSTRUM)
 	LDFLAGS = $(FFTW_LIBS) -mkl=sequential
 endif
 ifeq ($(MACHINE),SUPERMUC)
-	FC = h5pfc -fpp -r8 -O3
-	HDF5_LIBS = $(HDF5_F90_SHLIB) $(HDF5_SHLIB) -L$(SZIP_LIBDIR) -lz -ldl -lm
-	LDFLAGS = $(MKL_LIB) $(FFTW_LIB) $(HDF5_LIBS)
+	FC = mpif90 -fpp -r8 -O3 $(HDF5_INC)
+	LDFLAGS = $(FFTW_LIB) $(HDF5_F90_SHLIB) $(HDF5_SHLIB) -qmkl=sequential
 endif
 
 ifeq ($(MACHINE),SNELLIUS)
@@ -115,18 +110,23 @@ OBJS += obj/AddLatentHeat.o obj/DeallocatePFVariables.o obj/ExplicitTermsPhi.o \
 	obj/InterpTempMgrd.o obj/SolveImpEqnUpdate_Phi.o obj/CreateICPF.o \
 	obj/ImmersedBoundary.o
 
+# # Object files associated with the immersed boundary method
+OBJS += obj/SolveImpEqnUpdate_Temp_ibm.o obj/SolveImpEqnUpdate_X_ibm.o \
+	obj/SolveImpEqnUpdate_YZ_ibm.o obj/topogr_ibm.o obj/SolveImpEqnUpdate_Sal_ibm.o \
+	obj/DeallocateIBMVars.o
+
 # Object files for plane writing
 OBJS += obj/mean_zplane.o
 
 # Module object files
 MOBJS = obj/param.o obj/decomp_2d.o obj/AuxiliaryRoutines.o obj/decomp_2d_fft.o \
-	obj/HermiteInterpolations.o obj/GridModule.o obj/h5_tools.o obj/means.o
+	obj/HermiteInterpolations.o obj/GridModule.o obj/h5_tools.o obj/means.o obj/ibm_param.o
 
 #=======================================================================
 #  Files that create modules:
 #=======================================================================
 MFILES = param.F90 decomp_2d.F90 AuxiliaryRoutines.F90 decomp_2d_fft.F90 \
-	HermiteInterpolations.F90 GridModule.F90
+	HermiteInterpolations.F90 GridModule.F90 ibm_param.F90
 
 #============================================================================ 
 #  make PROGRAM   
@@ -148,7 +148,17 @@ $(OBJDIR)/AuxiliaryRoutines.o: src/flow_solver/AuxiliaryRoutines.F90
 $(OBJDIR)/decomp_2d.o: src/flow_solver/2decomp/decomp_2d.F90
 	$(FC) -c -o $@ $< $(LDFLAGS)
 $(OBJDIR)/decomp_2d_fft.o: src/flow_solver/2decomp/decomp_2d_fft.F90
-	$(FC) -c -o $@ $< $(LDFLAGS) 
+	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/ibm_param.o: src/ibm/ibm_param.F90
+	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/GridModule.o: src/flow_solver/GridModule.F90
+	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/HermiteInterpolations.o: src/multires/HermiteInterpolations.F90
+	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/h5_tools.o: src/h5tools/h5_tools.F90
+	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/means.o: src/h5tools/means.F90 obj/ibm_param.o
+	$(FC) -c -o $@ $< $(LDFLAGS)
 $(OBJDIR)/%.o: src/%.F90 $(MOBJS)
 	$(FC) -c -o $@ $< $(LDFLAGS)
 $(OBJDIR)/%.o: src/flow_solver/%.F90 $(MOBJS)
@@ -163,6 +173,8 @@ $(OBJDIR)/%.o: src/multires/phase-field/%.F90 $(MOBJS)
 	$(FC) -c -o $@ $< $(LDFLAGS)
 $(OBJDIR)/%.o: src/multires/salinity/%.F90 $(MOBJS)
 	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/%.o: src/ibm/%.F90 $(MOBJS)
+	$(FC) -c -o $@ $< $(LDFLAGS)
 
 #============================================================================
 #  Clean up 

docs/ibm.md

@@ -0,0 +1,102 @@
+# Immersed boundary method
+
+When using the immersed boundary method, we use the direct forcing method as described by [Fadlun et al. (2000)](https://doi.org/10.1006/jcph.2000.6484).
+This must be incorporated into the solution of the implicit step of the time stepper, and is taken care of by the subroutines `SolveImpEqnUpdateXX_ibm` stored in the `ibm` subdirectory.
+We distinguish points as either solid/fluid/boundary using the arrays `ibmaskXX`.
+
+## Modification to the implicit solver
+We show below the modifications to the implicit solve for the wall-normal velocity component $u$ (or `vx`) but the other components are analogous.
+In solving the implicit step, we solve the matrix system
+
+$$
+a_{ij} (\Delta u)_j = RHS^*_{i} ,
+$$
+
+where in the standard fluid domain $a_{ij}$ defines the wall-normal diffusion operator (see [the time stepper documentation](../numerics/#crank-nicolson-semi-implicit-diffusion) for more details).
+Recall that $\Delta u = u^{l+1} - u^l$ is the velocity increment over the RK substep.
+For all points away from the solid immersed boundaries, $RHS^*$ is simply the sum of the explicit terms and the explicit portion of the Crank-Nicolson term.
+
+### Solid points (`ibmask==0`)
+
+Inside the solid, we enforce $u^{l+1}=0$.
+If we consider the velocity at position $x_k$, this is achieved by setting $a_{kk}=1$, and $a_{kj}=0$ for $j\neq k$, with $RHS^*_k=-u^l_k$.
+
+$$
+(1)\Delta u = u^{l+1} - u^l = RHS^* = -u^l \quad \Rightarrow u^{l+1}=0 .
+$$
+
+Ideally, in the solid $u^l$ would be zero, but the pressure correction can introduce small velocities in the solid, so it is good practice to ensure the right hand side exactly cancels this previous velocity.
+
+### Boundary points (`ibmask==±1`)
+
+The boundary points are where we must take the most care.
+These points are defined as being those grid points adjacent (in the wall-normal direction `x`) to the immersed solids.
+For the points just above the interface, we set `ibmask=1` and for points just below a solid we set `ibmask=-1`, and we do not solve the Navier--Stokes equations.
+Rather, we interpolate the velocity linearly between the boundary and the second point into the liquid:
+
+$$
+u_k^{l+1} = \frac{\delta_2}{\delta_1 + \delta_2} u_{k+1}^{l+1} ,
+$$
+
+where $\delta_1$ is the $x$-distance between the first two grid points in the fluid, and $\delta_2$ is the $x$-distance between the first fluid grid point and the solid boundary.
+
+Since the interpolation involves the adjacent velocity at the next time step, this requires a further modification to the implicit solver.
+Rewriting in terms of the velocity increments, the condition we want to impose is
+
+$$
+\Delta u_k - \frac{\delta_2}{\delta_1 + \delta_2} \Delta u_{k+1} = \frac{\delta_2}{\delta_1 + \delta_2} u_{k+1}^l - u_k^l
+$$
+
+Note also that $k+1$ will be replaced by $k-1$ for boundary points where the liquid phase is below the solid boundary.
+The matrix for the implicit solver therefore reads
+
+$$
+a_{kj} = \begin{cases} 1 & j=k \\ -\delta_2/(\delta_1+\delta_2) & j= k+1 \\ 0 & \textrm{otherwise} \end{cases}
+$$
+
+with $RHS^*$ calculated as the right hand side of the equation above.
+This step is performed in the `SolveImpEqnUpdate*_ibm` routines, where the distance ratio in the above expressions is stored in the arrays `distx`, `disty`,... which are defined in `topogr_ibm`.
+
+### Temperature at the boundary
+
+The above treatment of the boundary points is a special case enforcing a fixed velocity of zero.
+For the temperature field, sometimes we want to enforce a specific non-zero temperature in the solid.
+This can be achieved simply by prescribing
+
+$$
+T_k = T_b + \frac{T_{k+1} - T_b}{\delta_1 + \delta_2} \delta_2 = \frac{\delta_2}{\delta_1 + \delta_2} T_{k+1} + \left( 1 - \frac{\delta_2}{\delta_1 + \delta_2} \right) T_b
+$$
+
+where $T_b$ is the fixed temperature enforced in the solid boundaries.
+Compared to the above case of zero velocity, all that is required is and addition of $T_b$ to the $RHS^*$ in the solid points, and an addition of $(1 - \delta_2/(\delta_1 + \delta_2)) T_b$ to the $RHS^*$ for the boundary points.
+
+### Concentration field: No-flux conditions and multi-resolution considerations
+
+Imposing zero-flux conditions at the immersed boundary is a bit more troublesome than enforcing a fixed value.
+Due to the complex nature of the three-dimensional semi-implicit time stepping on a staggered grid, we enforce zero gradients in *every* direction for points across the boundaries rather than only the normal gradient.
+In the wall-normal direction, we apply a similar approach to that above.
+
+A key difference here, is that the boundary points denoted `ibmask=±1` are now the closest points to the boundary *inside* the solid.
+Consider the point $x_k$ being one of these such boundary points.
+For the concentration field $C$, we then assume that $\partial_x C=0$ across the boundary by setting $C_{k}$ equal to $C_{k+1}$:
+
+$$
+C_k^{l+1} = C_{k+1}^{l+1} \ \Rightarrow \Delta C_k - \Delta C_{k+1} = C_{k+1}^l - C_k^l
+$$
+
+This is implemented in the same way as above into the implicit solver, just without any length ratios since they are not needed here.
+
+In the horizontal directions, there is no implicit step, so we cannot apply direct forcing for the immersed boundary method.
+To apply the zero-gradient condition, we locate the points adjacent to the boundary.
+At these points, we can write out the diffusive term as
+
+$$
+\kappa \partial_{yy} C = \frac{\kappa (C_{j+1} - C_j) - \kappa (C_j - C_{j-1})}{(\Delta y)^2}
+$$
+
+and set $C_{j-1}=C_j$ to eliminate half of the expression (where we consider point $y_{j-1}$ to be in the solid phase) and impose $\partial_y C=0$ across the boundary.
+For diffusion in the $x$-direction, we do not need to hard-code the modified diffusion since the direct immersed boundary forcing overrides the diffusion at the relevant points.
+
+Since the concentration field inside the solid has no physical meaning, we can do whatever is most convenient numerically for the points inside the solid (away from the boundary).
+This has no effect on the surrounding fluid region, since we have prescribed all the boundary points.
+In the current implementation, we choose to allow diffusion of $C$ inside the solid, such that the global solution is relatively smooth, which leads to a more reliable interpolation of $C$ onto the coarse velocity grid.

docs/numerics.md

@@ -63,19 +63,21 @@ $$
 \boldsymbol{\nabla}\cdot\boldsymbol{u}^{l+1} = 0 \approx \boldsymbol{\nabla}\cdot\boldsymbol{\hat{u}} - \alpha_l \Delta t \nabla^2 \Phi^{l+1} .
 $$
 
-The derivatives in the periodic directions are computed using Fourier transforms, whereas the wall-normal derivatives use a finite difference.
-Denoting $\boldsymbol{\tilde{u}}$ and $\tilde{\Phi}$ as the ($yz$-)Fourier transforms of $\boldsymbol{\hat{u}}$ and $\Phi^{l+1}$ respectively, we can write
+The divergence of $\hat{\boldsymbol{u}}$ is computed using finite differences, and then Fourier transformed.
+The derivatives for the pressure correction in the periodic directions can be expressed using Fourier transforms, whereas the wall-normal derivative must use a finite difference.
+Denoting $\boldsymbol{\tilde{u}}$ and $\tilde{\Phi}$ as the ($yz$-)Fourier transforms of $\boldsymbol{\hat{u}}$ and $\Phi^{l+1}$ respectively, we can now write the problem as a linear system
 
 $$
-\mathcal{D}\tilde{u} - ik_y \tilde{v} - ik_z\tilde{w} = \alpha_l \Delta t\left( \mathcal{DG}\tilde{\Phi} - (k_y^2 + k_z^2)\tilde{\Phi}\right) ,
+\alpha_l \Delta t\left[ \mathcal{DG} - (k_y^2 + k_z^2)\mathcal{I}\right]\tilde{\Phi} = \widetilde{\nabla\cdot \boldsymbol{u}},
 $$
 
 where $\mathcal{D}$ is the wall-normal discrete divergence operator.
-Note that since we are aiming to set the divergence of $\boldsymbol{u}^{l+1}$ to zero, we use the composite discrete operator $\mathcal{DG}$ for the second derivative of $\Phi$ rather than the discrete Laplacian $\mathcal{L}$ used above in the equations of motion.
-These are *not* the same operator.
-The above equation can be expressed as a tridiagonal matrix problem for each Fourier mode and solved using similar LAPACK routines as described above for the semi-implicit step.
+For each Fourier mode $k_y$, $k_z$, this provides a 1-D linear system that can be expressed as a tridiagonal matrix and solved for $\tilde{\Phi}$ using similar LAPACK routines as for the implicit step above.
 This is performed in the `SolvePressureCorrection` subroutine.
 
+
+Note that since we are aiming to set the divergence of $\boldsymbol{u}^{l+1}$ to zero, we use the composite discrete operator $\mathcal{DG}$ for the second derivative of $\Phi$ rather than the discrete Laplacian $\mathcal{L}$ used above in the equations of motion.
+These are *not* the same operator.
 For completeness, the composite operator $\mathcal{DG}$ takes the form
 
 $$

examples/HeatedWall/bou.in

@@ -52,4 +52,4 @@ xplusU  xminusU     dPdy    dPdz    MELT
 
 Phase-field parameters
 pf_D(*Pe)   pf_A    pf_S    pf_Tm   IBM     pf_IC
-0.1         10.0    10.0    0.0     1       1
+0.1         10.0    10.0    0.0     0       1

examples/PFWall/bou.in

@@ -52,4 +52,4 @@ xplusU  xminusU     dPdy    dPdz    MELT
 
 Phase-field parameters
 pf_D(*Pe)   pf_A    pf_S    pf_Tm   IBM     pf_IC
-0.48        1.0     10.0    0.0     1       1
+0.48        1.0     10.0    0.0     0       1

examples/RBPoiseuille/bou.in

@@ -52,4 +52,4 @@ xplusU  xminusU     dPdy    dPdz    MELT
 
 Phase-field parameters
 pf_A(*Pe)   pf_C    pf_S    pf_Tm   IBM     pf_IC
-0.1         10.0    10.0    0.0     1       1
+0.1         10.0    10.0    0.0     0       1

examples/RayleighBenardConvection/bou.in

@@ -52,4 +52,4 @@ xplusU  xminusU     dPdy    dPdz    MELT
 
 Phase-field parameters
 pf_D(*Pe)   pf_A    pf_S    pf_Tm   IBM     pf_IC
-0.1         10.0    10.0    0.0     1       1
+0.1         10.0    10.0    0.0     0       1

examples/StratifiedKH/bou.in

@@ -1,10 +1,10 @@
 Base grid resolution / number of substeps
 NXM     NYM     NZM     NSST(>1=RK,else AB)
-128     128     1       3
+192     256     1       3
 
 Refined grid: on/off and grid size
 MULTIRES    NXMR    NYMR    NZMR
-1           256     256     1
+1           384     512     1
 
 Salinity and Phase-field (on/off)
 FLAGSAL     FLAGPF
@@ -16,11 +16,11 @@ NREAD       IRESET
 
 Time steps / time limit (wall/flow) 
 NTST        WALLTIMEMAX     TMAX
-1000000     7000            500.0
+1000000     600            100.0
 
 Time interval for saving stats/movie frames
 TOUT    TFRAME  SAVE_3D
-0.1     0.1     50.0
+0.5     0.5     50.0
 
 Domain size (keep ALX3 = 1.0)
 ALX3   YLEN   ZLEN
@@ -52,4 +52,4 @@ xplusU  xminusU     dPdy    dPdz    MELT
 
 Phase-field parameters
 pf_D(*Pe)   pf_A    pf_S    pf_Tm   IBM     pf_IC
-0.1         10.0    10.0    0.0     1       1
+0.1         10.0    10.0    0.0     0       1

examples/ibm_test/bou.in

@@ -0,0 +1,55 @@
+Base grid resolution / number of substeps
+NXM     NYM     NZM     NSST(>1=RK,else AB)
+148     256     1       3
+
+Refined grid: on/off and grid size
+MULTIRES    NXMR    NYMR    NZMR
+1           592     1024     1
+
+Salinity and Phase-field (on/off)
+FLAGSAL     FLAGPF
+1           0
+
+Read flow / Reset logs
+NREAD       IRESET
+0           0
+
+Time steps / time limit (wall/flow) 
+NTST        WALLTIMEMAX     TMAX
+1000000     1200           100.0
+
+Time interval for saving stats/movie frames
+TOUT    TFRAME  SAVE_3D
+0.1     0.1    50.0
+
+Domain size (keep ALX3 = 1.0)
+ALX3   YLEN   ZLEN
+1.0    1.732    0.01
+
+Grid stretching parameters
+ISTR3   STR3    ISTR3R
+0       12.0    0
+
+Physical flow parameters / Free-fall timescale
+RAYT    PRAT    RAYS    PRAS    FFscaleS
+1e6     1.0     9.1e7   10.0    1
+
+Time stepping parameters
+IDTV    DT      RESID   CFLMAX  DTMIN   DTMAX
+1       1e-2    1.e-3   1.0     1e-8    1e-2
+
+Dirichlet/Neumann BC flag on upper(N) and lower(S) walls
+inslwS  inslwN  TfixS   TfixN   SfixS   SfixN
+1       1       1       1       0       0
+
+Active/passive scalar flags / Gravitational axis (x=1,y=2,z=3)
+active_T    active_S    gAxis
+0           1           1
+
+Wall velocity / Pressure gradient / Melt boundary condition
+xplusU  xminusU     dPdy    dPdz    MELT
+0.0     0.0         0.0     0.0     0
+
+Phase-field parameters
+pf_D(*Pe)   pf_A    pf_S    pf_Tm   IBM     pf_IC
+0.1         10.0    10.0    0.0     1       1

examples/mixed_convection/bou.in

@@ -52,4 +52,4 @@ xplusU  xminusU     dPdy    dPdz    MELT
 
 Phase-field parameters
 pf_A(*Pe)   pf_C    pf_S    pf_Tm   IBM     pf_IC
-0.1         10.0    10.0    0.0     1       1
+0.1         10.0    10.0    0.0     0       1