Skip to content

Allow strong pressure bcs to be used with stokes. #110

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Allow strong pressure bcs to be used with stokes. #110

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1134880
Allow strong pressure bcs to be used with stokes.
cianwilson Apr 21, 2016
eb2c367
Forgot to actually pass the inactive_mask through to the full project…
cianwilson Apr 21, 2016
74c0af5
Some improved comments.
cianwilson Apr 26, 2016
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
68 changes: 22 additions & 46 deletions assemble/Full_Projection.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ module Full_Projection

!--------------------------------------------------------------------------------------------------------------------
subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat,&
state, inner_mesh, auxiliary_matrix)
state, inner_mesh, auxiliary_matrix, inactive_mask)
!--------------------------------------------------------------------------------------------------------------------

! Solve Schur complement problem the nice way (using petsc) !!
Expand All @@ -84,6 +84,9 @@ subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat,&
type(mesh_type), intent(in):: inner_mesh
! p1-p1 stabilization matrix or free surface terms:
type(csr_matrix), optional, intent(in) :: auxiliary_matrix
! an array indicating which entries in the matrix are boundary conditions and
! therefore inactive (flagged as true)
logical, dimension(:), intent(in), optional :: inactive_mask

KSP ksp ! Object type for outer solve (i.e. A * delta_p = rhs)
Mat A ! PETSc Schur complement matrix (i.e. G^t*m^-1*G)
Expand All @@ -107,7 +110,7 @@ subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat,&
ewrite(2,*) 'Entering PETSc setup for Full Projection Solve'
call petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering,name,solver_option_path, &
lstartfromzero,inner_m,ctp_m,ct_m,x%option_path,pmat, &
rhs, state, inner_mesh, auxiliary_matrix)
rhs, state, inner_mesh, auxiliary_matrix, inactive_mask)

ewrite(2,*) 'Create RHS and solution Vectors in PETSc Format'
! create PETSc vec for rhs using above numbering:
Expand Down Expand Up @@ -140,7 +143,7 @@ end subroutine petsc_solve_full_projection
!--------------------------------------------------------------------------------------------------------
subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,solver_option_path, &
lstartfromzero,inner_m,div_matrix_comp, div_matrix_incomp,option_path,preconditioner_matrix,rhs, &
state, inner_mesh, auxiliary_matrix)
state, inner_mesh, auxiliary_matrix, inactive_mask)

!--------------------------------------------------------------------------------------------------------

Expand Down Expand Up @@ -181,6 +184,9 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so
! state, and inner_mesh are used to setup mg preconditioner of inner solve
type(state_type), intent(in):: state
type(mesh_type), intent(in):: inner_mesh
! an array indicating which entries in the matrix are boundary conditions and
! therefore inactive (flagged as true)
logical, dimension(:), intent(in), optional :: inactive_mask

type(vector_field), pointer :: positions, mesh_positions
type(petsc_csr_matrix), pointer:: rotation_matrix
Expand Down Expand Up @@ -209,7 +215,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so
integer reference_node, stat, i, rotation_stat
logical parallel, have_auxiliary_matrix, have_preconditioner_matrix

logical :: apply_reference_node, apply_reference_node_from_coordinates, reference_node_owned
integer :: j

! Sort option paths etc...
solver_option_path=complete_solver_option_path(option_path)
Expand All @@ -225,49 +231,18 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so
name=option_path
end if

! Are we applying a reference pressure node?
apply_reference_node = have_option(trim(option_path)//&
'/prognostic/reference_node')
apply_reference_node_from_coordinates = have_option(trim(option_path)//&
'/prognostic/reference_coordinates')

! If so, impose reference pressure node:
if(apply_reference_node) then

call get_option(trim(option_path)//&
'/prognostic/reference_node', reference_node)
if (GetProcNo()==1) then
ewrite(2,*) 'Imposing_reference_pressure_node'
allocate(ghost_nodes(1:1))
ghost_nodes(1) = reference_node
call set(rhs,reference_node,0.0) ! Modify RHS accordingly
else
allocate(ghost_nodes(1:0))
end if

elseif(apply_reference_node_from_coordinates) then

ewrite(1,*) 'Imposing_reference_pressure_node from user-specified coordinates'
positions => extract_vector_field(state, "Coordinate")
call find_reference_node_from_coordinates(positions,rhs%mesh,option_path,reference_node,reference_node_owned)
if(IsParallel()) then
if (reference_node_owned) then
allocate(ghost_nodes(1:1))
ghost_nodes(1) = reference_node
call set(rhs,reference_node,0.0)
else
allocate(ghost_nodes(1:0))
end if
else
allocate(ghost_nodes(1:1))
ghost_nodes(1) = reference_node
call set(rhs,reference_node,0.0)
end if

! create list of inactive, ghost_nodes
if(present(inactive_mask)) then
allocate( ghost_nodes(1:count(inactive_mask)) )
j=0
do i=1, size(inactive_mask)
if (inactive_mask(i)) then
j=j+1
ghost_nodes(j)=i
end if
end do
else

allocate(ghost_nodes(1:0))

allocate( ghost_nodes(1:0) )
end if

! Is auxiliary matrix present?
Expand All @@ -281,6 +256,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so
call allocate(petsc_numbering_p, &
nnodes=block_size(div_matrix_comp,1), nfields=1, &
halo=preconditioner_matrix%sparsity%row_halo, ghost_nodes=ghost_nodes)
deallocate(ghost_nodes)

! - why is this using the row halo of the preconditioner matrix when there might be rows missing?
! - same question about the nnodes use of the rows of the block of the divergence matrix?
Expand Down
16 changes: 14 additions & 2 deletions assemble/Momentum_Equation.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1806,6 +1806,10 @@ subroutine correct_pressure(state, prognostic_p_istate, x, u, p, old_p, delta_p,

type(vector_field), pointer :: positions

! An array indicating the location of the pressure reference node and/or
! boundary conditions. Only used for schur_complement solves.
logical, dimension(:), allocatable :: inactive_mask

ewrite(1,*) 'Entering correct_pressure'


Expand All @@ -1831,14 +1835,22 @@ subroutine correct_pressure(state, prognostic_p_istate, x, u, p, old_p, delta_p,

! Solve for the change in pressure, delta_p
if(full_schur) then
allocate(inactive_mask(size(ct_m(prognostic_p_istate)%ptr, 1)))
inactive_mask = .false.
call apply_dirichlet_conditions(inactive_mask, projec_rhs, p, dt=1.0/(dt*theta_pg*theta_divergence))
call impose_reference_pressure_node(inactive_mask, projec_rhs, positions, trim(p%option_path))

if(assemble_schur_auxiliary_matrix) then
call petsc_solve_full_projection(delta_p, ctp_m(prognostic_p_istate)%ptr, inner_m(prognostic_p_istate)%ptr, ct_m(prognostic_p_istate)%ptr, projec_rhs, &
full_projection_preconditioner, state(prognostic_p_istate), u%mesh, &
auxiliary_matrix=schur_auxiliary_matrix)
auxiliary_matrix=schur_auxiliary_matrix, inactive_mask=inactive_mask)
else
call petsc_solve_full_projection(delta_p, ctp_m(prognostic_p_istate)%ptr, inner_m(prognostic_p_istate)%ptr, ct_m(prognostic_p_istate)%ptr, projec_rhs, &
full_projection_preconditioner, state(prognostic_p_istate), u%mesh)
full_projection_preconditioner, state(prognostic_p_istate), u%mesh, &
inactive_mask=inactive_mask)
end if

deallocate(inactive_mask)
else
call petsc_solve(delta_p, cmc_m, projec_rhs, state(prognostic_p_istate))
end if
Expand Down
94 changes: 93 additions & 1 deletion femtools/Boundary_Conditions.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ module boundary_conditions
end interface

interface set_reference_node
module procedure set_reference_node_scalar, set_reference_node_vector_petsc
module procedure set_reference_node_scalar, set_reference_node_scalar_inactive_mask, set_reference_node_vector_petsc
end interface set_reference_node

interface has_boundary_condition
Expand All @@ -122,6 +122,7 @@ module boundary_conditions

interface apply_dirichlet_conditions
module procedure apply_dirichlet_conditions_scalar, &
apply_dirichlet_conditions_scalar_inactive_mask, &
apply_dirichlet_conditions_scalar_lumped, &
apply_dirichlet_conditions_vector, &
apply_dirichlet_conditions_vector_petsc_csr, &
Expand Down Expand Up @@ -1567,6 +1568,50 @@ subroutine set_reference_node_scalar(matrix, node, rhs, reference_value, referen

end subroutine set_reference_node_scalar

subroutine set_reference_node_scalar_inactive_mask(inactive_mask, node, rhs, reference_value, reference_node_owned)
!!< Sets a reference node for which the value is fixed in the equation
!!< This is typically done for a Poisson equation with all Neumann
!!< bcs to eliminate the spurious freedom of adding a constant value
!!< to the solution.
!!< This version of this routine returns the reference node location through
!!< an inactive mask array (inactive_mask), that is true in the location of
!!< the reference node.
logical, dimension(:), intent(inout) :: inactive_mask
integer, intent(in):: node
!! if rhs is not provided, you have to make sure the rhs at
!! the reference node has the right value, usually 0, yourself:
type(scalar_field), optional, intent(inout) :: rhs
!! by default the field gets set to 0 at the reference node
real, optional, intent(in) :: reference_value
!! in parallel all processes need to call this routine, only one
!! of them actually sets the reference node - this processor should
!! call with reference_node_owned=.true., other processes with reference_node_owned=.false.
!! if reference_node_owned is not present, we will assume that only process with rank 0
!! owns and thus sets the reference node
logical, optional, intent(in) :: reference_node_owned

logical:: lnode_owned

if (present(reference_node_owned)) then
lnode_owned = reference_node_owned
else
lnode_owned = GetProcNo()==1
end if

if (.not. lnode_owned) return

inactive_mask(node) = .true.

if (present(rhs)) then
if (present(reference_value)) then
call set(rhs, node, reference_value)
else
call set(rhs, node, 0.0)
end if
end if

end subroutine set_reference_node_scalar_inactive_mask

subroutine set_reference_node_vector_petsc(matrix, node, rhs, mask, reference_value, reference_node_owned)
!!< Sets a reference node for which the value is fixed in the equation
!!< This is typically done for a Poisson equation with all Neumann
Expand Down Expand Up @@ -1948,6 +1993,53 @@ subroutine apply_dirichlet_conditions_scalar(matrix, rhs, field, dt)

end subroutine apply_dirichlet_conditions_scalar

subroutine apply_dirichlet_conditions_scalar_inactive_mask(inactive_mask, rhs, field, dt)
!!< Apply dirichlet boundary conditions from field to the problem
!!<
!!< If dt is supplied, this assumes that boundary
!!< conditions are applied in rate of change form.
!!<
!!< This version of this routine returns the boundary condition location(s) through
!!< an inactive mask array (inactive_mask), that is true in the location of the
!!< dofs where the boundary condition is being applied.
logical, dimension(:):: inactive_mask
type(scalar_field), intent(inout) :: rhs
type(scalar_field), intent(in) :: field
real, optional, intent(in) :: dt

type(scalar_field), pointer:: surface_field
integer, dimension(:), pointer:: surface_node_list
character(len=FIELD_NAME_LEN):: bctype
integer :: i,j

bcloop: do i=1, get_boundary_condition_count(field)
call get_boundary_condition(field, i, type=bctype, &
surface_node_list=surface_node_list)

if (bctype/="dirichlet") cycle bcloop

inactive_mask(surface_node_list)=.true.

surface_field => extract_surface_field(field, i, "value")

if (present(dt)) then
do j=1, size(surface_node_list)
call set(rhs, surface_node_list(j), &
(node_val(surface_field, j)- &
node_val(field, surface_node_list(j)) &
)/dt)
end do
else
do j=1, size(surface_node_list)
call set(rhs, surface_node_list(j), &
node_val(surface_field, j))
end do
end if

end do bcloop

end subroutine apply_dirichlet_conditions_scalar_inactive_mask

subroutine apply_dirichlet_conditions_scalar_lumped(lhs, rhs, field, dt)
!!< Apply dirichlet boundary conditions from field to the problem
!!< defined by lhs and rhs. lhs is the diagonal of the normal matrix
Expand Down
58 changes: 56 additions & 2 deletions preprocessor/Boundary_Conditions_From_Options.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,10 @@ end subroutine projections_spherical_cartesian

end interface

interface impose_reference_pressure_node
module procedure impose_reference_pressure_node_matrix, impose_reference_pressure_node_inactive_mask
end interface impose_reference_pressure_node

contains

subroutine populate_boundary_conditions(states, suppress_warnings)
Expand Down Expand Up @@ -2342,7 +2346,7 @@ subroutine insert_flux_turbine_boundary_condition(field, bc_name, turbine_type)
end subroutine insert_flux_turbine_boundary_condition
end subroutine populate_flux_turbine_boundary_conditions

subroutine impose_reference_pressure_node(cmc_m, rhs, positions, option_path)
subroutine impose_reference_pressure_node_matrix(cmc_m, rhs, positions, option_path)
!!< If there are only Neumann boundaries on P, it is necessary to pin
!!< the value of the pressure at one point. As the rhs of the equation
!!< needs to be zeroed for this node, you will have to call this for
Expand Down Expand Up @@ -2387,7 +2391,57 @@ subroutine impose_reference_pressure_node(cmc_m, rhs, positions, option_path)

end if

end subroutine impose_reference_pressure_node
end subroutine impose_reference_pressure_node_matrix

subroutine impose_reference_pressure_node_inactive_mask(inactive_mask, rhs, positions, option_path)
!!< If there are only Neumann boundaries on P, it is necessary to pin
!!< the value of the pressure at one point. As the rhs of the equation
!!< needs to be zeroed for this node, you will have to call this for
!!< both pressure equations.
!!< This version of this routine returns the reference node location through
!!< an inactive mask array (inactive_mask), that is true in the location of
!!< the reference node.
logical, dimension(:), intent(inout) :: inactive_mask
type(scalar_field), intent(inout):: rhs
type(vector_field), intent(inout) :: positions
character(len=*), intent(in) :: option_path

integer :: stat, stat2
integer :: reference_node

logical :: apply_reference_node, apply_reference_node_from_coordinates, reference_node_owned

apply_reference_node = have_option(trim(complete_field_path(option_path, stat=stat2))//&
&"/reference_node")
apply_reference_node_from_coordinates = have_option(trim(complete_field_path(option_path, stat=stat2))//&
&"/reference_coordinates")

if(apply_reference_node) then

call get_option(trim(complete_field_path(option_path, stat=stat2))//&
&"/reference_node", reference_node, stat=stat)

if (stat==0) then
! all processors now have to call this routine, although only
! process 1 sets it
ewrite(1,*) 'Imposing_reference_pressure_node at node',reference_node
call set_reference_node(inactive_mask, reference_node, rhs)
end if

elseif(apply_reference_node_from_coordinates) then

ewrite(1,*) 'Imposing_reference_pressure_node at user-specified coordinates'
call find_reference_node_from_coordinates(positions,rhs%mesh,option_path,reference_node,reference_node_owned)

if(IsParallel()) then
call set_reference_node(inactive_mask, reference_node, rhs, reference_node_owned=reference_node_owned)
else
call set_reference_node(inactive_mask, reference_node, rhs)
end if

end if

end subroutine impose_reference_pressure_node_inactive_mask

subroutine find_reference_node_from_coordinates(positions,mesh,option_path,reference_node,reference_node_owned)
!! This routine determines which element contains the reference coordinates and,
Expand Down