Skip to main content

STDBE511 example 5

  • Both fluid and porous medium are considered
  • L = 0.1 m
  • dx = 0.1/50
  • K0 = E-10 m2
  • Da = E-8
  • V_TOP = 0.01
  • dt = dx / V_TOP = 0.2
  • StabOption=1

Mesh:

[TOC]

!!! note "Some simulation parameter"

#define _K0_ 1.0E-10
#define _PHI_ 0.5_DFP
#define _DT_ 0.2
#define _VTOP_ 0.01_DFP
#define _MU_FLUID_ 0.001
#define _RHO_FLUID_ 1000.0
#define _TIMESTEPS_ 100 ! maxI4B
#define _MAX_ITER_ 20
#define _tDBC4P_ 0
#define _REF_PRESSURE_NODE_ 90
#define _STAB_PARAM_OPTION_ 1
#define _MATERTIAL_INTERFACES_ [7]
#define _SAVE_STATE_FREQ_ 100
#define _WRITE_FREQ_ 1
varNameComment
_K0_permeability of the soil
_PHI_porosity of the soil
_DT_time step size
_VTOP_velocity at top
_MU_FLUID_dynamic viscosity of fluid
_RHO_FLUID_mass-density of fluid
_TIMESTEPS_number of time steps
_MAX_ITER_maximum number of Newton iteration
_tDBC4P_total Dirichlet BC for pressure
_REF_PRESSURE_NODE_reference node number for pressure
_STAB_PARAM_OPTION_stabilization option
_MATERTIAL_INTERFACES_mesh-id of fluid-soil interfaces
_SAVE_STATE_FREQ_save-state frequency
_WRITE_FREQ_write-frequency

!!! note "Include modules"

PROGRAM main
USE easifemBase
USE easifemClasses
USE easifemMaterials
USE easifemKernels
USE STDBE511_Class
IMPLICIT NONE

Declare variables

  TYPE(STDBE511_) :: obj
TYPE(ParameterList_) :: param
TYPE(HDF5File_) :: domainFile
TYPE( HDF5File_ ) :: outfile
TYPE( HDF5File_ ) :: statefile
TYPE(MeshSelection_) :: region
TYPE(Domain_), TARGET :: dom
CLASS(DirichletBC_), POINTER :: dbc => NULL()
CHARACTER(LEN=*), PARAMETER :: &
& domainFileName = "./mesh.h5"
CHARACTER( LEN = * ), PARAMETER :: &
& outfilePrefix="./output"
CHARACTER( LEN = * ), PARAMETER :: &
& statePrefix="./state"
INTEGER( I4B ), PARAMETER :: CoordinateSystem = KERNEL_2D
INTEGER( I4B ), PARAMETER :: SaveStateFrequency = _SAVE_STATE_FREQ_
INTEGER( I4B ), PARAMETER :: writeFrequency = _WRITE_FREQ_
INTEGER( I4B ) :: backupNumber
var-namecomment
objInstance of [[STDBE511_]] kernel
paramparameters instance of [[ParameterList_]]
domainFiledomain file [[HDF5File_]]
outfileoutput file [[HDF5File_]]
statefilestate file [[HDF5File_]]
regionmesh region [[MeshSelection_]]
domdomain [[Domain_]]
dbcinstance of dirichlet boundary condition [[NeumannBC_]]
domainFileNamedomain file name
outfilePrefixoutput file prefix
statePrefixstate file prefix
coordinateSystemcoordinate system
saveStateFrequencysaving frequency for state files
writeFrequencywrite frequency for output files
backupNumberinteral variable

!!! note "Linear solver parameters"

  INTEGER(I4B), PARAMETER :: LinSolver_name = LIS_GMRES
INTEGER(I4B), PARAMETER :: KrylovSubspaceSize = 50
INTEGER(I4B), PARAMETER :: maxIter_LinSolver = -1
REAL(DFP), PARAMETER :: rtol_LinSolver=1.0D-6
REAL(DFP), PARAMETER :: atol_LinSolver=1.0D-10
INTEGER(I4B), PARAMETER :: preconditionOption = NO_PRECONDITION
INTEGER(I4B), PARAMETER :: convergenceIn_LinSolver = convergenceInRes
INTEGER(I4B), PARAMETER :: convergenceType_LinSolver = relativeConvergence
LOGICAL( LGT ), PARAMETER :: relativeToRHS = .FALSE.
var-namecomment
LinSolver_namename of linear solver
KrylovSubspaceSizea parameter for linsol
maxIter_LinSolvermax number of iterations in linsol
rtol_LinSolverrelative tolerance for linsol
atol_LinSolverabsolute tolerance for linsol
preconditionOptionprecondition option
convergenceIn_LinSolverconvergence criteria for linsol
convergenceType_LinSolverconvergence criteria for linsol: relative or absolute.
relativeToRHSif convergence is relative then relative to RHS or not.

Currently, NO_PRECONDITON means diagonal preconditioning will be selected.

!!! note "Iteration parameters"

  INTEGER(I4B), PARAMETER :: nnt = 2
INTEGER( I4B ), PARAMETER :: stabParamOption = _STAB_PARAM_OPTION_
INTEGER( I4B ), PARAMETER :: maxIter= _MAX_ITER_
INTEGER(I4B), PARAMETER :: TotalTimeSteps = _TIMESTEPS_
REAL(DFP), PARAMETER :: startTime=0.0_DFP
REAL(DFP), PARAMETER :: dt= _DT_
REAL(DFP), PARAMETER :: endTime=dt*TotalTimeSteps
REAL( DFP ), PARAMETER :: atol_Pressure=1.0D-10
REAL( DFP ), PARAMETER :: rtol_Pressure=1.0D-10
REAL( DFP ), PARAMETER :: atol_Velocity=1.0D-6
REAL( DFP ), PARAMETER :: rtol_Velocity=1.0D-6
REAL( DFP ), PARAMETER :: tol_steadyState = 1.0E-8
LOGICAL( LGT ), PARAMETER :: resetIteration = .FALSE.
LOGICAL( LGT ), PARAMETER :: resetTimeStep = .FALSE.
var-namecomment
nntnumber of nodes in time
stabParamOptionstabilization parameter option
maxItermaximum number of iteration
TotalTimeStepstotal number of time steps for the simulation
startTimestart time of simulation
dttime step size of simulation
endTimeend time of simulation
atol_Pressureabsolute tolerance for checking convergence in pressure
rtol_Pressurerelative tolerance for checking convergence in pressure
atol_Velocityabsolute tolerance for checking convergence in velocity
rtol_Velocityrelative tolerance for checking convergence in velocity
tol_steadyStatetolerance for checking steady state
resetIterationreset variables to zero after each iteration
resetTimeStepreset variables to zero after each time step

!!! note "Materials and boundary conditions"

  INTEGER(I4B), PARAMETER :: tPorousMaterials = 2
INTEGER(I4B), PARAMETER :: tFluidMaterials = 1
INTEGER(I4B), PARAMETER :: tDirichletBCForVelocity = 3
INTEGER(I4B), PARAMETER :: tDirichletBCForPressure = _tDBC4P_
INTEGER( I4B ), PARAMETER :: refPressureNode = _REF_PRESSURE_NODE_
INTEGER( I4B ), PARAMETER :: materialInterfaces(1) = _MATERTIAL_INTERFACES_
REAL(DFP), PARAMETER :: porosity= _PHI_
REAL(DFP), PARAMETER :: porosity_fluid=1.0_DFP
REAL(DFP), PARAMETER :: permeability = _K0_
REAL(DFP), PARAMETER :: permeability_fluid=1.0E+30
REAL(DFP), PARAMETER :: massdensity_soil=1700.0
REAL(DFP), PARAMETER :: massdensity_fluid = _RHO_FLUID_
REAL(DFP), PARAMETER :: dynamicViscosity = _MU_FLUID_
REAL(DFP), PARAMETER :: V_TOP= _VTOP_
var-namecomment
tPorousMaterialstotal porous media
tFluidMaterialstotal fluid material
tDirichletBCForVelocitytotal dirichlet boundary condition for velocity
tDirichletBCForPressuretotal dirichlet boundary condition for pressure
refPressureNodereference node number for pressure
porosityporosity of soil
massdensity_soilmass density of soil
permeabilitypermeability of soil
massdensity_fluidmass density of fluid
dynamicViscositydynamic viscosity of fluid
V_TOPvelocity at the top

!!! note "Internal variables"

  LOGICAL( LGT ) :: convg
INTEGER( I4B ) :: iter
INTEGER( I4B ) :: its
REAL( DFP ) :: t1, t2, t2t1

SetSTDBE511Param

  • [[STDBE511_#SetSTDBE511Param]]
  CALL FPL_INIT(); CALL param%Initiate()
!!
CALL SetSTDBE511Param( &
& param=param, &
& materialInterfaces=materialInterfaces, &
& stabParamOption=stabParamOption, &
& engine="NATIVE_SERIAL", &
& nnt=nnt, &
& startTime=startTime, &
& endTime=endTime, &
& dt=dt, &
& coordinateSystem=coordinateSystem, &
& maxIter=maxIter, &
& atoleranceForPressure=atol_pressure, &
& atoleranceForVelocity=atol_velocity, &
& rtoleranceForPressure=rtol_pressure, &
& rtoleranceForVelocity=rtol_velocity, &
& toleranceForSteadyState=tol_steadyState, &
& tPorousMaterials=tPorousMaterials, &
& tFluidMaterials=tFluidMaterials, &
& tDirichletBCForVelocity=tDirichletBCForVelocity, &
& tDirichletBCForPressure=tDirichletBCForPressure, &
& domainFile=domainFileName)

SetLinSolverParam

  CALL SetLinSolverParam( &
& param=param, &
& solverName=LinSolver_name,&
& preconditionOption=preconditionOption, &
& convergenceIn=convergenceIn_LinSolver, &
& convergenceType=convergenceType_LinSolver, &
& maxIter=maxIter_LinSolver, &
& relativeToRHS=relativeToRHS, &
& KrylovSubspaceSize=KrylovSubspaceSize, &
& rtol=rtol_LinSolver, &
& atol=atol_LinSolver )

Domain Initiate

Domain for pressure and velocity field

  CALL domainFile%Initiate( filename=domainFileName, MODE="READ")
CALL domainFile%Open()
CALL dom%Initiate(domainFile, "")
CALL domainFile%Deallocate()

Initiate

Initiate the STDBE511 kernel [[STDBE511_#Initiate]]

  CALL obj%Initiate(param=param, dom=dom)

AddPorousMaterial

Add Porous Material (bottom region)

  • [[MeshSelection_#Initiate]]
  • [[MeshSelection_#Add]]
  • [[MeshSelection_#SetPorousMaterialParam]]
  CALL region%Initiate(isSelectionByMeshID=.TRUE.)
CALL region%Add(dim=obj%nsd, meshID=[1])
CALL region%Add(dim=obj%nsd-1, meshID=materialInterfaces )
  CALL SetPorousMaterialParam( &
& param=param, &
& name="porousMaterial", &
& massdensity=massdensity_soil, &
& porosity=porosity, &
& permeability=permeability )
  • [[STDBE511_#AddPorousMaterial]]
  CALL obj%AddPorousMaterial( &
& materialNo=1, &
& materialName="porousMaterial", &
& param=param, &
& region=region)

Deallocate the region, so that we can use it again.

  CALL region%Deallocate()

Adding another porous medium (top part)

Initiate region:

  CALL region%Initiate(isSelectionByMeshID=.TRUE.)
CALL region%Add(dim=obj%nsd, meshID=[2])

Set porous material properties.

  CALL SetPorousMaterialParam( &
& param=param, &
& name="porousMaterial", &
& massdensity=massdensity_fluid, &
& porosity=porosity_fluid, &
& permeability=permeability_fluid )

Add porous material to kernel, and deallocate the region.

  CALL obj%AddPorousMaterial( &
& materialNo=2, &
& materialName="porousMaterial", &
& param=param, &
& region=region)
!!
CALL region%Deallocate()

AddFluidMaterial

  • [[STDBE511_#AddFluidMaterial]]
  CALL region%Initiate(isSelectionByMeshID=.TRUE.)
CALL region%Add(dim=obj%nsd, meshID=[1,2])
CALL region%Add(dim=obj%nsd-1, meshID=materialInterfaces )
!!
CALL SetFluidMaterialParam( &
& param=param, &
& name="fluidMaterial", &
& massDensity=massdensity_fluid, &
& dynamicViscosity=dynamicViscosity )
!!
CALL obj%AddFluidMaterial( &
& materialNo=1, &
& materialName="fluidMaterial", &
& param=param, &
& region=region )
!!
CALL region%Deallocate()
!!

DirichletBoundaryCondition

In this problem we only have dirichlet boundary condition for velocity. At all boundary we have dirichlet boundary condition.

Vx=0

Select the mesh boundary where vx=0v_{x}=0 is prescribed.

  • [[MeshSelection_#Initiate]]
  • [[MeshSelection_#Add]]
  • [[MeshSelection_#Set]]
  CALL region%Initiate( &
& isSelectionByMeshID=.TRUE., &
& isSelectionByNodeNum=.TRUE. )
CALL region%Add(dim=obj%nsd-1, meshID=[1,2,6])
CALL region%Add( &
& nodenum=dom%getInternalNptrs( dim=obj%nsd-1, entityNum=[3,5] ) )
CALL region%Set()

Set the parameters for dirichlet boundary condition. [[NeumannBC_#SetDirichletBCParam]]

  CALL SetDirichletBCParam( &
& param=param, &
& name="Vx=0", &
& idof=1, &
& nodalValueType=Constant, &
& useFunction=.FALSE.)

Add dirichlet boundary condition to kernel [[STDBE511_#AddVelocityDirichletBC]], [[STDBE511_#GetVelocityDirichletBCPointer]]

  CALL obj%AddVelocityDirichletBC(dbcNo=1, param=param, boundary=region)
dbc => obj%GetVelocityDirichletBCPointer(dbcNo=1)
CALL dbc%Set(ConstantNodalValue=0.0_DFP)
dbc => NULL()
CALL region%Deallocate()

Vy=0

  CALL SetDirichletBCParam( &
& param=param, &
& name="Vy=0", &
& idof=2, &
& nodalValueType=Constant, &
& useFunction=.FALSE.)
  CALL region%Initiate(isSelectionByMeshID=.TRUE.)
CALL region%Add(dim=obj%nsd-1, meshID=[1,2,3,4,5,6])
CALL region%Set()
  CALL obj%AddVelocityDirichletBC(dbcNo=2, param=param, boundary=region)
dbc => obj%GetVelocityDirichletBCPointer(dbcNo=2)
CALL dbc%Set(ConstantNodalValue=0.0_DFP)
dbc => NULL()
CALL region%Deallocate()

Vx=Vtop

  CALL SetDirichletBCParam( &
& param=param, &
& name="Vx=V", &
& idof=1, &
& nodalValueType=Constant, &
& useFunction=.FALSE.)
  CALL region%Initiate(isSelectionByMeshID=.TRUE.)
CALL region%Add(dim=obj%nsd-1, meshID=[4])
CALL region%Set()
  CALL obj%AddVelocityDirichletBC(dbcNo=3, param=param, boundary=region)
dbc => obj%GetVelocityDirichletBCPointer(dbcNo=3)
CALL dbc%Set(ConstantNodalValue=V_TOP)
dbc => NULL()
CALL region%Deallocate()

Pressure boundary condition

Usually, we do not apply Dirichlet boundary condition for pressure.

  IF( tDirichletBCForPressure .GT. 0 ) THEN
CALL SetDirichletBCParam( &
& param=param, &
& name="P=0", &
& idof=1, &
& nodalValueType=Constant, &
& useFunction=.FALSE.)
CALL region%Initiate( &
& isSelectionByNodeNum=.TRUE. )
CALL region%Add( nodenum=[refPressureNode])
CALL region%Set()
!!
CALL obj%AddPressureDirichletBC(dbcNo=1, param=param, boundary=region)
CALL region%Deallocate()
dbc => obj%GetPressureDirichletBCPointer(dbcNo=1)
CALL dbc%Set(ConstantNodalValue=0.0_DFP)
dbc => NULL()
END IF

Set

CALL obj%Set( )

MainComputation

Related to output and backup

  backupNumber = 0
CALL outfile%Initiate(outfilePrefix//tostring(backupNumber)//'.h5', "NEW")
CALL outfile%Open()
CALL Display( "main debug-1" )

Loop-Start: TimeStep

  DO its = 1, TotalTimeSteps
CALL Display( 'TimeStep = ' // tostring(its) )
CALL obj%setCurrentTimeStep(its)

Loop-Start: Iteration

    DO iter = 1, maxIter
CALL Display( " Iteration = " // tostring(iter) )
CALL obj%setIterationNumber(iter)
CALL obj%ApplyDirichletBC()
CALL obj%AssembleTanmat()
CALL obj%AssembleRHS()
CALL obj%Assemble()
CALL obj%Solve()
convg = obj%isConverged()
CALL obj%UpdateIteration(reset=resetIteration)
IF( convg ) EXIT
END DO

Loop-End: Iteration

    IF( .not. convg ) THEN
CALL Display( "program stopped, no convergence" )
STOP
END IF
    IF( obj%isSteadyState() .OR. its .EQ. TotalTimeSteps ) THEN
CALL Display( 'STEADY STATE REACHED' )
CALL statefile%Initiate(statePrefix//tostring(its)//".h5", "NEW")
CALL statefile%Open()
CALL obj%WriteData(statefile, "")
CALL statefile%Deallocate()
EXIT
END IF
    CALL obj%Update(reset=resetTimeStep)
    IF( MOD( its, writeFrequency ) .EQ. 0 ) THEN
CALL obj%WriteData(outfile, "/"//tostring(its))
END IF
    IF( MOD( its, SaveStateFrequency ) .EQ. 0 ) THEN
backupNumber=backupNumber+1
CALL statefile%Initiate(statePrefix//tostring(its)//".h5", "NEW")
CALL statefile%Open()
CALL obj%WriteData(statefile, "")
CALL statefile%Deallocate()
!!
CALL outfile%Deallocate()
CALL outfile%Initiate( &
& outfilePrefix//tostring(backupNumber)//'.h5', "NEW")
CALL outfile%Open()
END IF

Loop-End: TimeStep

  END DO

!!! note "cleanup"

  CALL outfile%Deallocate()
CALL obj%Deallocate()
CALL dom%Deallocate()
CALL param%Deallocate(); CALL FPL_FINALIZE()
END PROGRAM main

Results