From 113488005a3cfb1912e118453f6a7965c3269925 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 21 Apr 2016 12:33:43 +0100 Subject: [PATCH 1/3] Allow strong pressure bcs to be used with stokes. This is done by providing interfaces to bc application onto an inactive mask, this is then passed into full projection and used to define the ghost nodes. To make this all happen in one place I've also done this for reference nodes. Needs tests. --- assemble/Full_Projection.F90 | 64 ++++---------- assemble/Momentum_Equation.F90 | 8 ++ femtools/Boundary_Conditions.F90 | 87 ++++++++++++++++++- .../Boundary_Conditions_From_Options.F90 | 55 +++++++++++- 4 files changed, 165 insertions(+), 49 deletions(-) diff --git a/assemble/Full_Projection.F90 b/assemble/Full_Projection.F90 index fc698da57c..cd3b5d9df9 100644 --- a/assemble/Full_Projection.F90 +++ b/assemble/Full_Projection.F90 @@ -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) !! @@ -84,6 +84,7 @@ 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 + 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) @@ -107,7 +108,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: @@ -140,7 +141,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) !-------------------------------------------------------------------------------------------------------- @@ -181,6 +182,7 @@ 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 + logical, dimension(:), intent(in), optional :: inactive_mask type(vector_field), pointer :: positions, mesh_positions type(petsc_csr_matrix), pointer:: rotation_matrix @@ -209,7 +211,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) @@ -225,49 +227,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? @@ -281,6 +252,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? diff --git a/assemble/Momentum_Equation.F90 b/assemble/Momentum_Equation.F90 index 59b20c8fd1..cf234de50d 100644 --- a/assemble/Momentum_Equation.F90 +++ b/assemble/Momentum_Equation.F90 @@ -1806,6 +1806,8 @@ subroutine correct_pressure(state, prognostic_p_istate, x, u, p, old_p, delta_p, type(vector_field), pointer :: positions + logical, dimension(:), allocatable :: inactive_mask + ewrite(1,*) 'Entering correct_pressure' @@ -1831,6 +1833,10 @@ 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))) + 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, & @@ -1839,6 +1845,8 @@ subroutine correct_pressure(state, prognostic_p_istate, x, u, p, old_p, delta_p, 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) end if + + deallocate(inactive_mask) else call petsc_solve(delta_p, cmc_m, projec_rhs, state(prognostic_p_istate)) end if diff --git a/femtools/Boundary_Conditions.F90 b/femtools/Boundary_Conditions.F90 index ef9ad1a606..15a84b2e79 100644 --- a/femtools/Boundary_Conditions.F90 +++ b/femtools/Boundary_Conditions.F90 @@ -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 @@ -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, & @@ -1567,6 +1568,47 @@ 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. + 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 @@ -1948,6 +1990,49 @@ 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. + 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 diff --git a/preprocessor/Boundary_Conditions_From_Options.F90 b/preprocessor/Boundary_Conditions_From_Options.F90 index ae6f1fc851..b8f7c1805b 100644 --- a/preprocessor/Boundary_Conditions_From_Options.F90 +++ b/preprocessor/Boundary_Conditions_From_Options.F90 @@ -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) @@ -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 @@ -2387,7 +2391,54 @@ 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. + 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, From eb2c3671bfa30965d06e407a831cbdb278f234f2 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 21 Apr 2016 16:19:54 +0100 Subject: [PATCH 2/3] Forgot to actually pass the inactive_mask through to the full projection. --- assemble/Momentum_Equation.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/assemble/Momentum_Equation.F90 b/assemble/Momentum_Equation.F90 index cf234de50d..df3f207798 100644 --- a/assemble/Momentum_Equation.F90 +++ b/assemble/Momentum_Equation.F90 @@ -1834,16 +1834,18 @@ 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) From 74c0af5bdc952552675efe90fdcc1a3b96c04062 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Tue, 26 Apr 2016 14:38:25 +0100 Subject: [PATCH 3/3] Some improved comments. --- assemble/Full_Projection.F90 | 4 ++++ assemble/Momentum_Equation.F90 | 2 ++ femtools/Boundary_Conditions.F90 | 7 +++++++ preprocessor/Boundary_Conditions_From_Options.F90 | 3 +++ 4 files changed, 16 insertions(+) diff --git a/assemble/Full_Projection.F90 b/assemble/Full_Projection.F90 index cd3b5d9df9..64efc5a0c8 100644 --- a/assemble/Full_Projection.F90 +++ b/assemble/Full_Projection.F90 @@ -84,6 +84,8 @@ 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) @@ -182,6 +184,8 @@ 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 diff --git a/assemble/Momentum_Equation.F90 b/assemble/Momentum_Equation.F90 index df3f207798..dc0637c2fb 100644 --- a/assemble/Momentum_Equation.F90 +++ b/assemble/Momentum_Equation.F90 @@ -1806,6 +1806,8 @@ 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' diff --git a/femtools/Boundary_Conditions.F90 b/femtools/Boundary_Conditions.F90 index 15a84b2e79..6f34df10bf 100644 --- a/femtools/Boundary_Conditions.F90 +++ b/femtools/Boundary_Conditions.F90 @@ -1573,6 +1573,9 @@ subroutine set_reference_node_scalar_inactive_mask(inactive_mask, node, rhs, ref !!< 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 @@ -1995,6 +1998,10 @@ subroutine apply_dirichlet_conditions_scalar_inactive_mask(inactive_mask, rhs, f !!< !!< 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 diff --git a/preprocessor/Boundary_Conditions_From_Options.F90 b/preprocessor/Boundary_Conditions_From_Options.F90 index b8f7c1805b..abac22bd76 100644 --- a/preprocessor/Boundary_Conditions_From_Options.F90 +++ b/preprocessor/Boundary_Conditions_From_Options.F90 @@ -2398,6 +2398,9 @@ subroutine impose_reference_pressure_node_inactive_mask(inactive_mask, rhs, posi !!< 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