From 0a0d98cb5e414d2a152a1581835257c640fbf29f Mon Sep 17 00:00:00 2001 From: Chris Howland Date: Tue, 9 Jan 2024 10:49:03 +0100 Subject: [PATCH 01/16] remove salinity noise from DDC IC --- src/salinity.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/salinity.F90 b/src/salinity.F90 index ab82ad41..8a2d8285 100644 --- a/src/salinity.F90 +++ b/src/salinity.F90 @@ -155,7 +155,7 @@ subroutine CreateInitialSalinity end do end do end do - call AddSalinityNoise(amp=5e-3, localised=.false.) + ! call AddSalinityNoise(amp=5e-3, localised=.false.) !! Stratified shear layer setup else if ((RayS < 0) .and. (RayT < 0)) then From 6fa9e04979c94a30ed94960781eecea9259e5a55 Mon Sep 17 00:00:00 2001 From: Chris Howland Date: Wed, 10 Jan 2024 11:17:35 +0100 Subject: [PATCH 02/16] add first (broken) implicit solve for direct forcing --- src/flow_solver/ImplicitAndUpdateVX.F90 | 3 + src/flow_solver/ImplicitAndUpdateVY.F90 | 3 + src/flow_solver/ImplicitAndUpdateVZ.F90 | 3 + src/flow_solver/TimeMarcher.F90 | 2 +- src/phasefield.F90 | 248 ++++++++++++++++++++++-- 5 files changed, 241 insertions(+), 18 deletions(-) diff --git a/src/flow_solver/ImplicitAndUpdateVX.F90 b/src/flow_solver/ImplicitAndUpdateVX.F90 index f85ed39f..e261ccf1 100644 --- a/src/flow_solver/ImplicitAndUpdateVX.F90 +++ b/src/flow_solver/ImplicitAndUpdateVX.F90 @@ -16,6 +16,7 @@ subroutine ImplicitAndUpdateVX use local_arrays, only: vx,rhs,rux,qcap,pr use decomp_2d, only: xstart,xend use ibm_param + use afid_phasefield, only: SolveImpEqnUpdate_X_pf implicit none integer :: jc,kc integer :: km,kp,ic @@ -73,6 +74,8 @@ subroutine ImplicitAndUpdateVX if (IBM) then call SolveImpEqnUpdate_X_ibm + elseif (phasefield) then + call SolveImpEqnUpdate_X_pf else call SolveImpEqnUpdate_X end if diff --git a/src/flow_solver/ImplicitAndUpdateVY.F90 b/src/flow_solver/ImplicitAndUpdateVY.F90 index 049d6ee2..986138ee 100644 --- a/src/flow_solver/ImplicitAndUpdateVY.F90 +++ b/src/flow_solver/ImplicitAndUpdateVY.F90 @@ -16,6 +16,7 @@ subroutine ImplicitAndUpdateVY use local_arrays, only: vy,ruy,pr,rhs,dph use decomp_2d, only: xstart,xend use ibm_param + use afid_phasefield, only: SolveImpEqnUpdate_VY_pf implicit none integer :: kc,jmm,jc,ic integer :: kpp,kmm @@ -85,6 +86,8 @@ subroutine ImplicitAndUpdateVY if (IBM) then call SolveImpEqnUpdate_YZ_ibm(vy,rhs,ibmasky,disty) + elseif (phasefield) then + call SolveImpEqnUpdate_VY_pf else call SolveImpEqnUpdate_YZ(vy,rhs) end if diff --git a/src/flow_solver/ImplicitAndUpdateVZ.F90 b/src/flow_solver/ImplicitAndUpdateVZ.F90 index c2c05e7f..657ce76f 100644 --- a/src/flow_solver/ImplicitAndUpdateVZ.F90 +++ b/src/flow_solver/ImplicitAndUpdateVZ.F90 @@ -16,6 +16,7 @@ subroutine ImplicitAndUpdateVZ use local_arrays, only: vz,dq,ruz,rhs,pr use decomp_2d, only: xstart,xend use ibm_param + use afid_phasefield, only: SolveImpEqnUpdate_VZ_pf implicit none integer :: kc,jc,ic,imm integer :: kmm,kpp @@ -70,6 +71,8 @@ subroutine ImplicitAndUpdateVZ if (IBM) then call SolveImpEqnUpdate_YZ_ibm(vz,rhs,ibmaskz,distz) + else if (phasefield) then + call SolveImpEqnUpdate_VZ_pf else call SolveImpEqnUpdate_YZ(vz,rhs) end if diff --git a/src/flow_solver/TimeMarcher.F90 b/src/flow_solver/TimeMarcher.F90 index 5a5cc11e..3b3db22b 100644 --- a/src/flow_solver/TimeMarcher.F90 +++ b/src/flow_solver/TimeMarcher.F90 @@ -103,7 +103,7 @@ subroutine TimeMarcher ! iF(ANY(IsNaN(temp))) write(*,*)nrank,'NaN in TEMP post-implicit',ns ! if (phasefield .and. IBM) call ImmersedBoundary - if (phasefield .and. (pf_direct_force > 0)) call ForceIceVelZero + ! if (phasefield .and. (pf_direct_force > 0)) call ForceIceVelZero call update_halo(vy,lvlhalo) call update_halo(vz,lvlhalo) diff --git a/src/phasefield.F90 b/src/phasefield.F90 index b1390a82..44792893 100644 --- a/src/phasefield.F90 +++ b/src/phasefield.F90 @@ -98,8 +98,9 @@ subroutine CreateInitialPhase !! Favier (2019) appendix A.3 validation case elseif (pf_IC==3) then - call set_flat_interface(0.5, .true.) - call add_temperature_mode(amp=0.1, ymode=2, zmode=2) + h0 = 0.4 + call set_flat_interface(h0, .true.) + call add_temperature_mode(amp=1e-2, ymode=10, zmode=2, h0=h0) else if (pf_IC==4) then call set_ice_sphere(r0=0.1) @@ -237,34 +238,35 @@ end subroutine set_multicomponent_interface !> Add a modal perturbation to the temperature field in the lower half of the domain !! Intended for use with phase-field validation case of 2D Melting RBC from Favier et al (2019) -!! Vertical structure is of the form sin^2(2 pi x) -subroutine add_temperature_mode(amp, ymode, zmode) +!! Vertical structure is of the form sin^2(pi x/h0) +subroutine add_temperature_mode(amp, ymode, zmode, h0) use local_arrays, only: temp real, intent(in) :: amp !! amplitude of perturbation integer, intent(in) :: ymode !! mode number in y integer, intent(in) :: zmode !! mode number in z (only used if nzm>1) - integer :: i, j, k, kmid + real, intent(in) :: h0 !! initial interface height + integer :: i, j, k real :: xxx, yyy, zzz - kmid = nxm/2 do i=xstart(3),xend(3) + zzz = zm(i)/zlen do j=xstart(2),xend(2) + yyy = ym(j)/ylen if (nzm > 1) then - do k=1,kmid ! If domain 3D, add in z perturbation too + do k=1,nxm ! If domain 3D, add in z perturbation too xxx = xm(k) - yyy = ym(j)/ylen - zzz = zm(i)/zlen - - temp(k,j,i) = temp(k,j,i) & - + amp*sin(2.0*pi*ymode*yyy)*cos(2.0*pi*zmode*zzz)*sin(2.0*pi*xxx)**2 + if (xxx < h0) then + temp(k,j,i) = temp(k,j,i) & + + amp*sin(2.0*pi*ymode*yyy)*cos(2.0*pi*zmode*zzz)*sin(pi*xxx/h0)**2 + end if end do else - do k=1,kmid + do k=1,nxm xxx = xm(k) - yyy = ym(j)/ylen - - temp(k,j,i) = temp(k,j,i) & - + amp*sin(2.0*pi*ymode*yyy)*sin(2.0*pi*xxx)**2 + if (xxx < h0) then + temp(k,j,i) = temp(k,j,i) & + + amp*sin(2.0*pi*ymode*yyy)*sin(pi*xxx/h0)**2 + end if end do end if end do @@ -839,4 +841,216 @@ subroutine CreatePhaseH5Groups(filename) end subroutine CreatePhaseH5Groups +!> Implicit solve for vx when using direct forcing constraint +subroutine SolveImpEqnUpdate_X_pf + use local_arrays, only: vx, rhs + integer :: ic, jc, kc + real :: fkl(1:nx) ! RHS vector for implicit solve + real :: philoc ! local value of phi on vx-grid + real :: amkl(1:nxm) ! Lower diagonal of LHS matrix + real :: ackl(1:nx) ! Diagonal of LHS matrix + real :: apkl(1:nxm) ! Upper diagonal of LHS matrix + real :: appk(1:nxm-1) + integer :: ipkv(1:nx), info + real :: ackl_b, betadx + + betadx = beta*al + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + + amkl(:) = 0.0 + ackl(:) = 1.0 + apkl(:) = 0.0 + + fkl(1) = -vx(1,jc,ic) + + do kc=2,nxm + philoc = 0.5*(phic(kc,jc,ic) + phic(kc+1,jc,ic)) + if (philoc > pf_direct_force) then ! Solid phase + amkl(kc-1) = 0.d0 + ackl(kc) = 1.d0 + apkl(kc) = 0.d0 + fkl(kc) = -vx(kc,jc,ic) + else ! Liquid phase + ackl_b=1.0d0/(1.0d0-ac3ck(kc)*betadx) + amkl(kc-1)=-am3ck(kc)*betadx*ackl_b + ackl(kc)=1.d0 + apkl(kc)=-ap3ck(kc)*betadx*ackl_b + fkl(kc) = rhs(kc,jc,ic)*ackl_b + end if + end do + + fkl(nx) = -vx(nx,jc,ic) + + call dgttrf(nx, amkl, ackl, apkl, appk, ipkv, info) + call dgttrs('N',nx,1,amkl,ackl,apkl,appk,ipkv,fkl,nx,info) + + do kc=2,nxm + vx(kc,jc,ic) = vx(kc,jc,ic) + fkl(kc) + end do + + end do + end do + +end subroutine SolveImpEqnUpdate_X_pf + +!> Implicit solve for vy, setting zero in solid +subroutine SolveImpEqnUpdate_VY_pf + use local_arrays, only: vy, rhs + integer :: ic, jc, kc + real :: fkl(nxm) ! RHS vector for implicit solve + real :: philoc ! local value of phi on vx-grid + real :: amkl(nxm-1) ! Lower diagonal of LHS matrix + real :: ackl(nxm) ! Diagonal of LHS matrix + real :: apkl(nxm-1) ! Upper diagonal of LHS matrix + integer :: ipkv(1:nxm), info + real :: ackl_b, betadx + real :: appk(nx) + + betadx = beta*al + appk(:) = 0.0 + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + + amkl(:) = 0.0 + ackl(:) = 1.0 + apkl(:) = 0.0 + + do kc=1,nxm + philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc+1,ic)) + if (philoc > pf_direct_force) then ! Solid phase + amkl(kc-1) = 0.d0 + ackl(kc) = 1.d0 + apkl(kc) = 0.d0 + fkl(kc) = -vy(kc,jc,ic) + else ! Liquid phase + ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) + if (kc > 1) amkl(kc-1)=-am3sk(kc)*betadx*ackl_b + ackl(kc)=1.d0 + if (kc < nxm) apkl(kc)=-ap3sk(kc)*betadx*ackl_b + fkl(kc) = rhs(kc,jc,ic)*ackl_b + end if + end do + + call dgttrf(nxm, amkl, ackl, apkl, appk, ipkv, info) + call dgttrs('N',nxm,1,amkl,ackl,apkl,appk,ipkv,fkl,nxm,info) + + do kc=1,nxm + vy(kc,jc,ic) = vy(kc,jc,ic) + fkl(kc) + end do + end do + end do + +end subroutine SolveImpEqnUpdate_VY_pf + + +!> Implicit solve for vz, setting zero in solid +subroutine SolveImpEqnUpdate_VZ_pf + use local_arrays, only: vz, rhs + integer :: ic, jc, kc + real :: fkl(nxm) ! RHS vector for implicit solve + real :: philoc ! local value of phi on vx-grid + real :: amkl(nxm-1) ! Lower diagonal of LHS matrix + real :: ackl(nxm) ! Diagonal of LHS matrix + real :: apkl(nxm-1) ! Upper diagonal of LHS matrix + integer :: ipkv(1:nxm), info + real :: ackl_b, betadx + real :: appk(nx-3) + + betadx = beta*al + + do ic=xstart(3),xend(3) + do jc=xstart(2),xend(2) + + amkl(:) = 0.0 + ackl(:) = 1.0 + apkl(:) = 0.0 + + do kc=1,nxm + philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic+1)) + if (philoc > pf_direct_force) then ! Solid phase + amkl(kc-1) = 0.d0 + ackl(kc) = 1.d0 + apkl(kc) = 0.d0 + fkl(kc) = -vz(kc,jc,ic) + else ! Liquid phase + ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) + if (kc > 1) amkl(kc-1)=-am3sk(kc)*betadx*ackl_b + ackl(kc)=1.d0 + if (kc < nxm) apkl(kc)=-ap3sk(kc)*betadx*ackl_b + fkl(kc) = rhs(kc,jc,ic)*ackl_b + end if + end do + + call dgttrf(nxm, amkl, ackl, apkl, appk, ipkv, info) + call dgttrs('N',nxm,1,amkl,ackl,apkl,appk,ipkv,fkl,nxm,info) + + do kc=1,nxm + vz(kc,jc,ic) = vz(kc,jc,ic) + fkl(kc) + end do + end do + end do + +end subroutine SolveImpEqnUpdate_VZ_pf + +! !> Implicit solve for vy, vz, setting zero in solid +! subroutine SolveImpEqnUpdate_YZ_pf(q, rhs, axis) +! real, intent(inout) :: q(1:nx,xstart(2)-lvlhalo:xend(2)+lvlhalo, & +! & xstart(3)-lvlhalo:xend(3)+lvlhalo) !> variable to update +! real, intent(inout) :: rhs(1:nx,xstart(2):xend(2),xstart(3):xend(3)) !> rhs to update with +! character, intent(in) :: axis !> component axis of velocity (to identify grid stagger) +! integer :: ic, jc, kc +! real :: fkl(nxm) ! RHS vector for implicit solve +! real :: philoc ! local value of phi on vx-grid +! real :: amkl(nxm-1) ! Lower diagonal of LHS matrix +! real :: ackl(nxm) ! Diagonal of LHS matrix +! real :: apkl(nxm-1) ! Upper diagonal of LHS matrix +! integer :: ipkv(1:nxm), info +! real :: ackl_b, betadx +! real :: appl(nx-3) + +! betadx = beta*al + +! do ic=xstart(3),xend(3) +! do jc=xstart(2),xend(2) + +! amkl(:) = 0.0 +! ackl(:) = 1.0 +! apkl(:) = 0.0 + +! do kc=1,nxm +! if (axis=="y") then +! philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc+1,ic)) +! elseif (axis=="z") then +! philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic+1)) +! else +! philoc = phic(kc,jc,ic) +! end if +! if (philoc > pf_direct_force) then ! Solid phase +! amkl(kc-1) = 0.d0 +! ackl(kc) = 1.d0 +! apkl(kc) = 0.d0 +! fkl(kc) = -q(kc,jc,ic) +! else ! Liquid phase +! ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) +! if (kc > 1) amkl(kc-1)=-am3sk(kc)*betadx*ackl_b +! ackl(kc)=1.d0 +! if (kc < nxm) apkl(kc)=-ap3sk(kc)*betadx*ackl_b +! fkl(kc) = rhs(kc,jc,ic)*ackl_b +! end if +! end do + +! call dgttrf(nxm, amkl, ackl, apkl, appl, ipkv, info) +! call dgttrs('N',nxm,1,amkl,ackl,apkl,appl,ipkv,fkl,nxm,info) + +! do kc=1,nxm +! q(kc,jc,ic) = q(kc,jc,ic) + fkl(kc) +! end do +! end do +! end do + +! end subroutine SolveImpEqnUpdate_YZ_pf + end module afid_phasefield \ No newline at end of file From 6c49322c51f5b1533aba3c283f0b60554ab06d71 Mon Sep 17 00:00:00 2001 From: Chris Howland Date: Wed, 10 Jan 2024 11:18:39 +0100 Subject: [PATCH 03/16] update params for melt RBC test case --- examples/2DMeltingRBC/bou.in | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/2DMeltingRBC/bou.in b/examples/2DMeltingRBC/bou.in index 291eae50..b3a0e093 100644 --- a/examples/2DMeltingRBC/bou.in +++ b/examples/2DMeltingRBC/bou.in @@ -12,15 +12,15 @@ FLAGSAL FLAGPF Read flow / Reset logs NREAD IRESET -1 0 +0 0 Time steps / time limit (wall/flow) NTST WALLTIMEMAX TMAX -1000000 120 100.0 +1000000 3500 100.0 Time interval for saving stats/movie frames TOUT TFRAME SAVE_3D -0.01 0.01 10.0 +0.1 1.0 10.0 Domain size (keep ALX3 = 1.0) ALX3 YLEN ZLEN @@ -36,7 +36,7 @@ RAYT PRAT RAYS PRAS FFscaleS Time stepping parameters IDTV DT RESID CFLMAX DTMIN DTMAX -1 1e-3 1.e-3 1.0 1e-8 5e-3 +1 1e-3 1.e-3 1.0 1e-8 1e-3 Dirichlet/Neumann BC flag on upper(N) and lower(S) walls inslwS inslwN TfixS TfixN SfixS SfixN From bc7ccdab7f8bf0042f1cd62e361f34dbe58bf658 Mon Sep 17 00:00:00 2001 From: Chris Howland Date: Wed, 10 Jan 2024 13:55:32 +0100 Subject: [PATCH 04/16] fix index issue with implicit solve arrays --- src/phasefield.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/phasefield.F90 b/src/phasefield.F90 index 44792893..c3b668ca 100644 --- a/src/phasefield.F90 +++ b/src/phasefield.F90 @@ -921,9 +921,9 @@ subroutine SolveImpEqnUpdate_VY_pf do kc=1,nxm philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc+1,ic)) if (philoc > pf_direct_force) then ! Solid phase - amkl(kc-1) = 0.d0 + if (kc > 1) amkl(kc-1) = 0.d0 ackl(kc) = 1.d0 - apkl(kc) = 0.d0 + if (kc < nxm) apkl(kc) = 0.d0 fkl(kc) = -vy(kc,jc,ic) else ! Liquid phase ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) @@ -971,9 +971,9 @@ subroutine SolveImpEqnUpdate_VZ_pf do kc=1,nxm philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic+1)) if (philoc > pf_direct_force) then ! Solid phase - amkl(kc-1) = 0.d0 + if (kc > 1) amkl(kc-1) = 0.d0 ackl(kc) = 1.d0 - apkl(kc) = 0.d0 + if (kc < nxm) apkl(kc) = 0.d0 fkl(kc) = -vz(kc,jc,ic) else ! Liquid phase ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) @@ -1009,7 +1009,7 @@ end subroutine SolveImpEqnUpdate_VZ_pf ! real :: apkl(nxm-1) ! Upper diagonal of LHS matrix ! integer :: ipkv(1:nxm), info ! real :: ackl_b, betadx -! real :: appl(nx-3) +! real :: appk(nxm-2) ! betadx = beta*al @@ -1029,9 +1029,9 @@ end subroutine SolveImpEqnUpdate_VZ_pf ! philoc = phic(kc,jc,ic) ! end if ! if (philoc > pf_direct_force) then ! Solid phase -! amkl(kc-1) = 0.d0 +! if (kc > 1) amkl(kc-1) = 0.d0 ! ackl(kc) = 1.d0 -! apkl(kc) = 0.d0 +! if (kc < nxm) apkl(kc) = 0.d0 ! fkl(kc) = -q(kc,jc,ic) ! else ! Liquid phase ! ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) @@ -1042,8 +1042,8 @@ end subroutine SolveImpEqnUpdate_VZ_pf ! end if ! end do -! call dgttrf(nxm, amkl, ackl, apkl, appl, ipkv, info) -! call dgttrs('N',nxm,1,amkl,ackl,apkl,appl,ipkv,fkl,nxm,info) +! call dgttrf(nxm, amkl, ackl, apkl, appk, ipkv, info) +! call dgttrs('N',nxm,1,amkl,ackl,apkl,appk,ipkv,fkl,nxm,info) ! do kc=1,nxm ! q(kc,jc,ic) = q(kc,jc,ic) + fkl(kc) From 9658312a91871c4b365c44e35d6e082322c8b1fb Mon Sep 17 00:00:00 2001 From: Chris Howland Date: Fri, 26 Jan 2024 16:09:40 +0100 Subject: [PATCH 05/16] modify initial temp fields for PF --- src/flow_solver/CreateInitialConditions.F90 | 114 ++------------------ 1 file changed, 7 insertions(+), 107 deletions(-) diff --git a/src/flow_solver/CreateInitialConditions.F90 b/src/flow_solver/CreateInitialConditions.F90 index 14e2c294..c97751d7 100644 --- a/src/flow_solver/CreateInitialConditions.F90 +++ b/src/flow_solver/CreateInitialConditions.F90 @@ -14,7 +14,7 @@ subroutine CreateInitialConditions use decomp_2d, only: xstart,xend use mpih use afid_salinity, only: RayS - use afid_phasefield, only: pf_eps, read_phase_field_params + use afid_phasefield, only: pf_eps, read_phase_field_params, pf_Tm implicit none integer :: j,k,i,kmid real :: xxx,yyy,zzz,eps,varptb,amp @@ -250,86 +250,20 @@ subroutine CreateInitialConditions end if if (phasefield) then + ! Most of this is now in `afid_phasefield` in the routine `CreateInitialPhase` - if (pf_IC == 1) then ! 1D moving interface example - h0 = 0.1 ! Freezing if RAYT < 0 - Lambda = 0.620063 ! Melting if RAYT > 0 - t0 = pect * (h0/2/Lambda)**2 + if (pf_IC==3) then + h0 = 0.4 do i=xstart(3),xend(3) do j=xstart(2),xend(2) do k=1,nxm xxx = xm(k) + ! Piecewise linear base profile for Purseed et al if (xxx < h0) then - temp(k,j,i) = erf(xxx*sqrt(pect/t0)/2)/erf(Lambda) + temp(k,j,i) = 1.0 - (1.0 - pf_Tm)*xxx/h0 else - temp(k,j,i) = 1.0 - end if - if (RAYT > 0) temp(k,j,i) = 1.0 - temp(k,j,i) - end do - end do - end do - - else if (pf_IC == 2) then - do i=xstart(3),xend(3) - do j=xstart(2),xend(2) - do k=1,nxm - r = sqrt((xm(k) - 0.5)**2 + (ym(j) - ylen/2)**2) - temp(k,j,i) = 0.5*(1.0 + tanh(100.0*(r - 0.1))) - end do - end do - end do - - else if (pf_IC == 3) then ! Favier et al (2019) Appendix A3 Validation Case - eps = 0.1 - kmid = nxm/2 - do i=xstart(3),xend(3) - do j=xstart(2),xend(2) - if (nzm > 1) then - do k=1,kmid ! If domain 3D, add in z perturbation too - xxx = xm(k) - yyy = ym(j) - zzz = zm(i) - - temp(k,j,i) = temp(k,j,i) & - + eps*sin(4.0*pi*yyy)*cos(4.0*pi*zzz)*sin(2.0*pi*xxx)**2 - end do - else - do k=1,kmid - xxx = xm(k) - yyy = ym(j) - - temp(k,j,i) = temp(k,j,i) & - + eps*sin(4.0*pi*yyy)*sin(2.0*pi*xxx)**2 - end do - end if - end do - end do - - elseif (pf_IC == 4) then ! Ice block to compare with Neufeld et al (2010) - do i=xstart(3),xend(3) ! Ice at bottom of domain if RayT<0, at the top if RayT>0 - do j=xstart(2),xend(2) - r = sqrt((ym(j) - ylen/2.0)**2 + (zm(i) - zlen/2.0)**2) - do k=1,nxm - temp(k,j,i) = 1.0 - 0.25*(1.0 + sign(1.0,RayT)*tanh((xm(k) - alx3/2.0)/2.0/pf_eps)) & - *(1.0 - tanh((r - alx3/2.0)/2.0/pf_eps)) - end do - end do - end do - - else if (pf_IC == 5) then ! 1D supercooling example - h0 = 0.02 - Lambda = 0.060314 - t0 = pect * (h0/2/Lambda)**2 - do i=xstart(3),xend(3) - do j=xstart(2),xend(2) - do k=1,nxm - xxx = xm(k) - if (xxx > h0) then - temp(k,j,i) = erfc(xxx*sqrt(pect/t0)/2.0)/erfc(Lambda) - else - temp(k,j,i) = 1.0 + temp(k,j,i) = pf_Tm*(1.0 - xxx)/(1.0 - h0) end if - ! if (RAYT > 0) temp(k,j,i) = 1.0 - temp(k,j,i) end do end do end do @@ -416,21 +350,6 @@ subroutine CreateInitialConditions end if end if - ! if (IBM) then - ! do i=xstart(3),xend(3) - ! do j=xstart(2),xend(2) - ! h0 = 0.25 + (ym(j) - 0.5)**2 - ! do k=1,nxm - ! if (xm(k) < h0) then - ! temp(k,j,i) = 0.0 - ! else - ! temp(k,j,i) = 1.0 - ! end if - ! end do - ! end do - ! end do - ! end if - end if if (melt) then @@ -446,24 +365,5 @@ subroutine CreateInitialConditions end do end if - ! FAVIER ET AL. (2019) APPENDIX A1 VALIDATION CASE - ! if (phasefield) then - ! ! kmid = nxm/2 - ! eps = 8.041 - ! do i=xstart(3),xend(3) - ! do j=xstart(2),xend(2) - ! do k=1,nxm - ! temp(k,j,i) = (exp(-eps*(xm(k) - 1.0)) - 1.0)/(exp(eps) - 1.0) - ! end do - ! ! do k=1,kmid - ! ! temp(k,j,i) = 1.0 - 2.0*xm(k) - ! ! end do - ! ! do k=kmid+1,nxm - ! ! temp(k,j,i) = 0.1 - 0.2*xm(k) - ! ! end do - ! end do - ! end do - ! end if - return end subroutine CreateInitialConditions \ No newline at end of file From c3bcd7854c6efd9dc02c672963bc4c8e66d152e3 Mon Sep 17 00:00:00 2001 From: chowland Date: Tue, 13 Feb 2024 14:58:51 +0100 Subject: [PATCH 06/16] remove unnecessary linking flags --- Makefile | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/Makefile b/Makefile index 5a01f543..002b134c 100644 --- a/Makefile +++ b/Makefile @@ -144,49 +144,49 @@ $(PROGRAM): $(MOBJS) $(OBJS) # Dependencies #============================================================================ $(OBJDIR)/param.o: src/flow_solver/param.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/AuxiliaryRoutines.o: src/flow_solver/AuxiliaryRoutines.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/decomp_2d.o: src/flow_solver/2decomp/decomp_2d.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/decomp_2d_fft.o: src/flow_solver/2decomp/decomp_2d_fft.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/ibm_param.o: src/ibm/ibm_param.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/grid.o: src/grid.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/pressure.o: src/pressure.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/HermiteInterpolations.o: src/multires/HermiteInterpolations.F90 obj/ibm_param.o - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/h5_tools.o: src/h5tools/h5_tools.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/means.o: src/h5tools/means.F90 obj/ibm_param.o - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/IBMTools.o: src/ibm/IBMTools.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/salinity.o: src/salinity.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/phasefield.o: src/phasefield.F90 obj/salinity.o - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/moisture.o: src/moisture.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/time_averaging.o: src/time_averaging.F90 - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/spectra.o: src/spectra.F90 obj/time_averaging.o obj/pressure.o - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/%.o: src/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/%.o: src/flow_solver/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/%.o: src/h5tools/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/%.o: src/multires/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/%.o: src/multires/IC_interpolation/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< $(OBJDIR)/%.o: src/ibm/%.F90 $(MOBJS) - $(FC) -c -o $@ $< $(LDFLAGS) + $(FC) -c -o $@ $< #============================================================================ # Clean up From b85bad1c84178980a290d0271dea2075da422aaf Mon Sep 17 00:00:00 2001 From: chowland Date: Tue, 13 Feb 2024 15:00:16 +0100 Subject: [PATCH 07/16] update supermuc packages --- Makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Makefile b/Makefile index 002b134c..6c5ae89c 100644 --- a/Makefile +++ b/Makefile @@ -8,7 +8,7 @@ FLAVOUR=GNU # Intel: 2022 intel/2022a FFTW/3.3.10-GCC-11.3.0 HDF5/1.12.2-iimpi-2021a # IRENE (Intel): flavor/hdf5/parallel hdf5 fftw3/gnu # MARENOSTRUM (Intel): fabric intel mkl impi hdf5 fftw szip -# SUPERMUC (Intel): fftw hdf5 +# SUPERMUC (Intel): spack/23.1.0 intel-toolkit/2023.1.0 fftw hdf5 # DISCOVERER: # GNU: hdf5/1/1.14/latest-gcc-openmpi fftw/3/latest-gcc-openmpi lapack # Intel: hdf5/1/1.14/latest-intel-openmpi fftw/3/latest-gcc-openmpi mkl @@ -65,7 +65,7 @@ ifeq ($(MACHINE),MARENOSTRUM) LDFLAGS = $(FFTW_LIBS) -mkl=sequential endif ifeq ($(MACHINE),SUPERMUC) - FC = mpif90 -fpp -r8 -O3 $(HDF5_INC) + FC = mpiifort -r8 -O3 $(HDF5_INC) LDFLAGS = $(FFTW_LIB) $(HDF5_F90_SHLIB) $(HDF5_SHLIB) -qmkl=sequential endif From 87258a2f16bf3855aab75a68477f4b693a85afb0 Mon Sep 17 00:00:00 2001 From: chowland Date: Tue, 5 Mar 2024 17:33:42 +0100 Subject: [PATCH 08/16] use generic yz implicit solve for pf --- src/flow_solver/ImplicitAndUpdateVY.F90 | 4 +- src/flow_solver/ImplicitAndUpdateVZ.F90 | 4 +- src/phasefield.F90 | 132 +++--------------------- 3 files changed, 20 insertions(+), 120 deletions(-) diff --git a/src/flow_solver/ImplicitAndUpdateVY.F90 b/src/flow_solver/ImplicitAndUpdateVY.F90 index 986138ee..7e2e3a71 100644 --- a/src/flow_solver/ImplicitAndUpdateVY.F90 +++ b/src/flow_solver/ImplicitAndUpdateVY.F90 @@ -16,7 +16,7 @@ subroutine ImplicitAndUpdateVY use local_arrays, only: vy,ruy,pr,rhs,dph use decomp_2d, only: xstart,xend use ibm_param - use afid_phasefield, only: SolveImpEqnUpdate_VY_pf + use afid_phasefield, only: SolveImpEqnUpdate_YZ_pf implicit none integer :: kc,jmm,jc,ic integer :: kpp,kmm @@ -87,7 +87,7 @@ subroutine ImplicitAndUpdateVY if (IBM) then call SolveImpEqnUpdate_YZ_ibm(vy,rhs,ibmasky,disty) elseif (phasefield) then - call SolveImpEqnUpdate_VY_pf + call SolveImpEqnUpdate_YZ_pf(vy,rhs,'y') else call SolveImpEqnUpdate_YZ(vy,rhs) end if diff --git a/src/flow_solver/ImplicitAndUpdateVZ.F90 b/src/flow_solver/ImplicitAndUpdateVZ.F90 index 657ce76f..a1c36ae9 100644 --- a/src/flow_solver/ImplicitAndUpdateVZ.F90 +++ b/src/flow_solver/ImplicitAndUpdateVZ.F90 @@ -16,7 +16,7 @@ subroutine ImplicitAndUpdateVZ use local_arrays, only: vz,dq,ruz,rhs,pr use decomp_2d, only: xstart,xend use ibm_param - use afid_phasefield, only: SolveImpEqnUpdate_VZ_pf + use afid_phasefield, only: SolveImpEqnUpdate_YZ_pf implicit none integer :: kc,jc,ic,imm integer :: kmm,kpp @@ -72,7 +72,7 @@ subroutine ImplicitAndUpdateVZ if (IBM) then call SolveImpEqnUpdate_YZ_ibm(vz,rhs,ibmaskz,distz) else if (phasefield) then - call SolveImpEqnUpdate_VZ_pf + call SolveImpEqnUpdate_YZ_pf(vz,rhs,'z') else call SolveImpEqnUpdate_YZ(vz,rhs) end if diff --git a/src/phasefield.F90 b/src/phasefield.F90 index c3b668ca..b503aeba 100644 --- a/src/phasefield.F90 +++ b/src/phasefield.F90 @@ -895,9 +895,12 @@ subroutine SolveImpEqnUpdate_X_pf end subroutine SolveImpEqnUpdate_X_pf -!> Implicit solve for vy, setting zero in solid -subroutine SolveImpEqnUpdate_VY_pf - use local_arrays, only: vy, rhs +!> Implicit solve for vy, vz, setting zero in solid +subroutine SolveImpEqnUpdate_YZ_pf(q, rhs, axis) + real, intent(inout) :: q(1:nx,xstart(2)-lvlhalo:xend(2)+lvlhalo, & + & xstart(3)-lvlhalo:xend(3)+lvlhalo) !> variable to update + real, intent(inout) :: rhs(1:nx,xstart(2):xend(2),xstart(3):xend(3)) !> rhs to update with + character, intent(in) :: axis !> component axis of velocity (to identify grid stagger) integer :: ic, jc, kc real :: fkl(nxm) ! RHS vector for implicit solve real :: philoc ! local value of phi on vx-grid @@ -906,10 +909,9 @@ subroutine SolveImpEqnUpdate_VY_pf real :: apkl(nxm-1) ! Upper diagonal of LHS matrix integer :: ipkv(1:nxm), info real :: ackl_b, betadx - real :: appk(nx) + real :: appk(nxm-2) betadx = beta*al - appk(:) = 0.0 do ic=xstart(3),xend(3) do jc=xstart(2),xend(2) @@ -919,62 +921,18 @@ subroutine SolveImpEqnUpdate_VY_pf apkl(:) = 0.0 do kc=1,nxm - philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc+1,ic)) - if (philoc > pf_direct_force) then ! Solid phase - if (kc > 1) amkl(kc-1) = 0.d0 - ackl(kc) = 1.d0 - if (kc < nxm) apkl(kc) = 0.d0 - fkl(kc) = -vy(kc,jc,ic) - else ! Liquid phase - ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) - if (kc > 1) amkl(kc-1)=-am3sk(kc)*betadx*ackl_b - ackl(kc)=1.d0 - if (kc < nxm) apkl(kc)=-ap3sk(kc)*betadx*ackl_b - fkl(kc) = rhs(kc,jc,ic)*ackl_b + if (axis=="y") then + philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc+1,ic)) + elseif (axis=="z") then + philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic+1)) + else + philoc = phic(kc,jc,ic) end if - end do - - call dgttrf(nxm, amkl, ackl, apkl, appk, ipkv, info) - call dgttrs('N',nxm,1,amkl,ackl,apkl,appk,ipkv,fkl,nxm,info) - - do kc=1,nxm - vy(kc,jc,ic) = vy(kc,jc,ic) + fkl(kc) - end do - end do - end do - -end subroutine SolveImpEqnUpdate_VY_pf - - -!> Implicit solve for vz, setting zero in solid -subroutine SolveImpEqnUpdate_VZ_pf - use local_arrays, only: vz, rhs - integer :: ic, jc, kc - real :: fkl(nxm) ! RHS vector for implicit solve - real :: philoc ! local value of phi on vx-grid - real :: amkl(nxm-1) ! Lower diagonal of LHS matrix - real :: ackl(nxm) ! Diagonal of LHS matrix - real :: apkl(nxm-1) ! Upper diagonal of LHS matrix - integer :: ipkv(1:nxm), info - real :: ackl_b, betadx - real :: appk(nx-3) - - betadx = beta*al - - do ic=xstart(3),xend(3) - do jc=xstart(2),xend(2) - - amkl(:) = 0.0 - ackl(:) = 1.0 - apkl(:) = 0.0 - - do kc=1,nxm - philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic+1)) if (philoc > pf_direct_force) then ! Solid phase if (kc > 1) amkl(kc-1) = 0.d0 ackl(kc) = 1.d0 if (kc < nxm) apkl(kc) = 0.d0 - fkl(kc) = -vz(kc,jc,ic) + fkl(kc) = -q(kc,jc,ic) else ! Liquid phase ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) if (kc > 1) amkl(kc-1)=-am3sk(kc)*betadx*ackl_b @@ -988,69 +946,11 @@ subroutine SolveImpEqnUpdate_VZ_pf call dgttrs('N',nxm,1,amkl,ackl,apkl,appk,ipkv,fkl,nxm,info) do kc=1,nxm - vz(kc,jc,ic) = vz(kc,jc,ic) + fkl(kc) + q(kc,jc,ic) = q(kc,jc,ic) + fkl(kc) end do end do end do -end subroutine SolveImpEqnUpdate_VZ_pf - -! !> Implicit solve for vy, vz, setting zero in solid -! subroutine SolveImpEqnUpdate_YZ_pf(q, rhs, axis) -! real, intent(inout) :: q(1:nx,xstart(2)-lvlhalo:xend(2)+lvlhalo, & -! & xstart(3)-lvlhalo:xend(3)+lvlhalo) !> variable to update -! real, intent(inout) :: rhs(1:nx,xstart(2):xend(2),xstart(3):xend(3)) !> rhs to update with -! character, intent(in) :: axis !> component axis of velocity (to identify grid stagger) -! integer :: ic, jc, kc -! real :: fkl(nxm) ! RHS vector for implicit solve -! real :: philoc ! local value of phi on vx-grid -! real :: amkl(nxm-1) ! Lower diagonal of LHS matrix -! real :: ackl(nxm) ! Diagonal of LHS matrix -! real :: apkl(nxm-1) ! Upper diagonal of LHS matrix -! integer :: ipkv(1:nxm), info -! real :: ackl_b, betadx -! real :: appk(nxm-2) - -! betadx = beta*al - -! do ic=xstart(3),xend(3) -! do jc=xstart(2),xend(2) - -! amkl(:) = 0.0 -! ackl(:) = 1.0 -! apkl(:) = 0.0 - -! do kc=1,nxm -! if (axis=="y") then -! philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc+1,ic)) -! elseif (axis=="z") then -! philoc = 0.5*(phic(kc,jc,ic) + phic(kc,jc,ic+1)) -! else -! philoc = phic(kc,jc,ic) -! end if -! if (philoc > pf_direct_force) then ! Solid phase -! if (kc > 1) amkl(kc-1) = 0.d0 -! ackl(kc) = 1.d0 -! if (kc < nxm) apkl(kc) = 0.d0 -! fkl(kc) = -q(kc,jc,ic) -! else ! Liquid phase -! ackl_b=1.0d0/(1.0d0-ac3sk(kc)*betadx) -! if (kc > 1) amkl(kc-1)=-am3sk(kc)*betadx*ackl_b -! ackl(kc)=1.d0 -! if (kc < nxm) apkl(kc)=-ap3sk(kc)*betadx*ackl_b -! fkl(kc) = rhs(kc,jc,ic)*ackl_b -! end if -! end do - -! call dgttrf(nxm, amkl, ackl, apkl, appk, ipkv, info) -! call dgttrs('N',nxm,1,amkl,ackl,apkl,appk,ipkv,fkl,nxm,info) - -! do kc=1,nxm -! q(kc,jc,ic) = q(kc,jc,ic) + fkl(kc) -! end do -! end do -! end do - -! end subroutine SolveImpEqnUpdate_YZ_pf +end subroutine SolveImpEqnUpdate_YZ_pf end module afid_phasefield \ No newline at end of file From 358a6fc08513f1b2f50b12e6065059362c52ba73 Mon Sep 17 00:00:00 2001 From: chowland Date: Tue, 5 Mar 2024 17:34:01 +0100 Subject: [PATCH 09/16] clean up old commands --- src/flow_solver/TimeMarcher.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/flow_solver/TimeMarcher.F90 b/src/flow_solver/TimeMarcher.F90 index 3b3db22b..65bdc7cf 100644 --- a/src/flow_solver/TimeMarcher.F90 +++ b/src/flow_solver/TimeMarcher.F90 @@ -86,7 +86,6 @@ subroutine TimeMarcher ! iF(ANY(IsNaN(temp))) write(*,*)nrank,'NaN in TEMP pre-implicit',ns if (phasefield) call update_halo(phi,lvlhalo) - ! if (phasefield .and. IBM) call UpdateIBMLocation call ImplicitAndUpdateVX call ImplicitAndUpdateVY @@ -102,9 +101,6 @@ subroutine TimeMarcher ! iF(ANY(IsNaN(vz))) write(*,*)nrank,'NaN in VZ post-implicit',ns ! iF(ANY(IsNaN(temp))) write(*,*)nrank,'NaN in TEMP post-implicit',ns - ! if (phasefield .and. IBM) call ImmersedBoundary - ! if (phasefield .and. (pf_direct_force > 0)) call ForceIceVelZero - call update_halo(vy,lvlhalo) call update_halo(vz,lvlhalo) From f77cfb7f04630a476f22a49fc3cd8a5c1956061b Mon Sep 17 00:00:00 2001 From: chowland Date: Wed, 6 Mar 2024 10:24:24 +0100 Subject: [PATCH 10/16] correct IC for 1D melting case --- src/phasefield.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/phasefield.F90 b/src/phasefield.F90 index b503aeba..d3df0bf4 100644 --- a/src/phasefield.F90 +++ b/src/phasefield.F90 @@ -89,7 +89,11 @@ subroutine CreateInitialPhase !! (RayT > 0: melting; RayT < 0: freezing) if (pf_IC==1) then h0 = 0.1 - call set_flat_interface(h0, .false.) + if (RayT > 0) then + call set_flat_interface(h0, .true.) + else + call set_flat_interface(h0, .false.) + end if call set_temperature_interface(h0, .false.) !! 2D axisymmetric melting example (disc radius 0.1) From 16e43543c9e2e2f90588878d6d75c98df5a029ea Mon Sep 17 00:00:00 2001 From: chowland Date: Wed, 6 Mar 2024 11:56:50 +0100 Subject: [PATCH 11/16] update axisym melt time step --- examples/2DAxisymMelting/bou.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/2DAxisymMelting/bou.in b/examples/2DAxisymMelting/bou.in index 19a1e199..aacbd248 100644 --- a/examples/2DAxisymMelting/bou.in +++ b/examples/2DAxisymMelting/bou.in @@ -36,7 +36,7 @@ RAYT PRAT RAYS PRAS FFscaleS Time stepping parameters IDTV DT RESID CFLMAX DTMIN DTMAX -1 1e-3 1.e-3 1.0 1e-8 1e-2 +0 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 From eab30e24d3b28957607e9c673438e86c14a316ea Mon Sep 17 00:00:00 2001 From: chowland Date: Wed, 6 Mar 2024 12:01:32 +0100 Subject: [PATCH 12/16] remove unused examples --- examples/2DIceBlock/bou.in | 55 -------------------------------------- examples/3DIceBlock/bou.in | 55 -------------------------------------- 2 files changed, 110 deletions(-) delete mode 100644 examples/2DIceBlock/bou.in delete mode 100644 examples/3DIceBlock/bou.in diff --git a/examples/2DIceBlock/bou.in b/examples/2DIceBlock/bou.in deleted file mode 100644 index 8f34bd47..00000000 --- a/examples/2DIceBlock/bou.in +++ /dev/null @@ -1,55 +0,0 @@ -Base grid resolution / number of substeps -NXM NYM NZM NSST(>1=RK,else AB) -128 256 1 3 - -Refined grid: on/off and grid size -MULTIRES NXMR NYMR NZMR -1 256 512 1 - -Salinity and Phase-field (on/off) -FLAGSAL FLAGPF -0 1 - -Read flow / Reset logs -NREAD IRESET -0 1 - -Time steps / time limit (wall/flow) -NTST WALLTIMEMAX TMAX -1000000 1200 50.0 - -Time interval for saving stats/movie frames -TOUT TFRAME SAVE_3D -0.1 0.1 10.0 - -Domain size (keep ALX3 = 1.0) -ALX3 YLEN ZLEN -1.0 2.0 0.01 - -Grid stretching parameters -ISTR3 STR3 ISTR3R -0 12.0 0 - -Physical flow parameters / Free-fall timescale -RAYT PRAT RAYS PRAS FFscaleS -1e6 10.0 1e8 1.0 0 - -Time stepping parameters -IDTV DT RESID CFLMAX DTMIN DTMAX -1 1e-3 1.e-3 1.0 1e-8 2e-3 - -Dirichlet/Neumann BC flag on upper(N) and lower(S) walls -inslwS inslwN TfixS TfixN SfixS SfixN -1 1 1 0 1 1 - -Active/passive scalar flags / Gravitational axis (x=1,y=2,z=3) -active_T active_S gAxis -1 0 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 -1.2e-1 10.0 1.0 0.0 0 4 \ No newline at end of file diff --git a/examples/3DIceBlock/bou.in b/examples/3DIceBlock/bou.in deleted file mode 100644 index 1893e36e..00000000 --- a/examples/3DIceBlock/bou.in +++ /dev/null @@ -1,55 +0,0 @@ -Base grid resolution / number of substeps -NXM NYM NZM NSST(>1=RK,else AB) -128 256 256 3 - -Refined grid: on/off and grid size -MULTIRES NXMR NYMR NZMR -1 256 512 512 - -Salinity and Phase-field (on/off) -FLAGSAL FLAGPF -0 1 - -Read flow / Reset logs -NREAD IRESET -0 1 - -Time steps / time limit (wall/flow) -NTST WALLTIMEMAX TMAX -1000000 1200 50.0 - -Time interval for saving stats/movie frames -TOUT TFRAME SAVE_3D -0.1 0.1 10.0 - -Domain size (keep ALX3 = 1.0) -ALX3 YLEN ZLEN -1.0 2.0 2.0 - -Grid stretching parameters -ISTR3 STR3 ISTR3R -0 12.0 0 - -Physical flow parameters / Free-fall timescale -RAYT PRAT RAYS PRAS FFscaleS --1e6 10.0 1e8 1.0 0 - -Time stepping parameters -IDTV DT RESID CFLMAX DTMIN DTMAX -1 1e-3 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 0 1 1 1 - -Active/passive scalar flags / Gravitational axis (x=1,y=2,z=3) -active_T active_S gAxis -1 0 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 -1.2e-1 10.0 1.0 0.0 0 4 \ No newline at end of file From 63cd86d6d5df8633719bf7d282ed38da0f5c7ad6 Mon Sep 17 00:00:00 2001 From: chowland Date: Wed, 6 Mar 2024 12:02:24 +0100 Subject: [PATCH 13/16] fix typo --- docs/inputs.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/inputs.md b/docs/inputs.md index 15f81e9f..2fd11094 100644 --- a/docs/inputs.md +++ b/docs/inputs.md @@ -115,4 +115,4 @@ The input parameter `ibm` can be used to specify different shapes of solid objec - `ibm = 4` produces a hexagonal lattice of spheres Of these options, only 1 and 4 have been reliably tested. -Perform your own validation casses before using. \ No newline at end of file +Perform your own validation cases before using. \ No newline at end of file From ffa554add58513295ec509cee006aaef8cd89ccd Mon Sep 17 00:00:00 2001 From: chowland Date: Wed, 6 Mar 2024 15:03:14 +0100 Subject: [PATCH 14/16] remove old forcing subroutine --- src/phasefield.F90 | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/src/phasefield.F90 b/src/phasefield.F90 index d3df0bf4..9f0e3154 100644 --- a/src/phasefield.F90 +++ b/src/phasefield.F90 @@ -380,34 +380,6 @@ subroutine ExplicitPhase end subroutine ExplicitPhase -!> Force the velocity field to be zero in the ice phase -!! Apply this only above the phase-field level pf_direct_force -subroutine ForceIceVelZero - use local_arrays, only: vx, vy, vz - integer :: i, j, k, im, jm - - do i=xstart(3),xend(3) - im = i - 1 - do j=xstart(2),xend(2) - jm = j - 1 - do k=2,nxm - if (0.5*(phic(k,j,i) + phic(k-1,j,i)) > pf_direct_force) then - vx(k,j,i) = 0.0 - end if - end do - do k=1,nxm - if (0.5*(phic(k,j,i) + phic(k,jm,i)) > pf_direct_force) then - vy(k,j,i) = 0.0 - end if - if (0.5*(phic(k,j,i) + phic(k,j,im)) > pf_direct_force) then - vz(k,j,i) = 0.0 - end if - end do - end do - end do - -end subroutine ForceIceVelZero - !> Compute the implicit terms for the phase-field evolution subroutine ImplicitPhase integer :: jc,kc,ic From f1c48598e82c0c0fe34abfda1600973585f3f428 Mon Sep 17 00:00:00 2001 From: chowland Date: Wed, 6 Mar 2024 15:24:59 +0100 Subject: [PATCH 15/16] fix math expression --- docs/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/index.md b/docs/index.md index 94457c70..924a0b40 100644 --- a/docs/index.md +++ b/docs/index.md @@ -21,7 +21,7 @@ Differences with [van der Poel et al (2015)](https://doi.org/10.1016/j.compfluid Differences with [Ostilla-Monico et al (2015)](https://doi.org/10.1016/j.jcp.2015.08.031): - Since the current code is modified from AFiD, it is pencil-parallelized in the periodic directions. This contrasts with the previous multi-resolution code, which was slab-parallelized in the wall-normal direction. -- For reasons linked to this change in parallelization, the only terms calculated implicitly are the diffusive terms with derivatives in the wall-normal direction (e.g. $`\nu \partial_{xx}u`$). All other terms are computed explicitly. +- For reasons linked to this change in parallelization, the only terms calculated implicitly are the diffusive terms with derivatives in the wall-normal direction (e.g. $\nu \partial_{xx}u$). All other terms are computed explicitly. - The multiple resolution strategy in time from [Ostilla-Monico et al (2015)](https://doi.org/10.1016/j.jcp.2015.08.031) is not yet implemented in the code. For now, we only rely on a CFL condition. ## Prerequisites From b1d3ea9fa62fe4f7d18d893fce861014e86fb124 Mon Sep 17 00:00:00 2001 From: chowland Date: Wed, 6 Mar 2024 15:25:11 +0100 Subject: [PATCH 16/16] switch recommendation to openmpi --- docs/prerequisites.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/prerequisites.md b/docs/prerequisites.md index a21cf1f1..b6aa5706 100644 --- a/docs/prerequisites.md +++ b/docs/prerequisites.md @@ -17,7 +17,7 @@ These instructions also work for machines running Ubuntu in the Windows Subsyste ## Single-line prerequisite installation ``` -sudo apt install build-essential gfortran libblas-dev liblapack-dev libhdf5-mpich-dev libfftw3-dev +sudo apt install build-essential gfortran libblas-dev liblapack-dev libhdf5-openmpi-dev libfftw3-dev ``` ## Using a new HPC?