diff --git a/assemble/Full_Projection.F90 b/assemble/Full_Projection.F90 index fc698da57c..64efc5a0c8 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,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) @@ -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: @@ -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) !-------------------------------------------------------------------------------------------------------- @@ -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 @@ -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) @@ -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? @@ -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? diff --git a/assemble/Momentum_Equation.F90 b/assemble/Momentum_Equation.F90 index 59b20c8fd1..dc0637c2fb 100644 --- a/assemble/Momentum_Equation.F90 +++ b/assemble/Momentum_Equation.F90 @@ -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' @@ -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 diff --git a/femtools/Boundary_Conditions.F90 b/femtools/Boundary_Conditions.F90 index ef9ad1a606..6f34df10bf 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,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 @@ -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 diff --git a/preprocessor/Boundary_Conditions_From_Options.F90 b/preprocessor/Boundary_Conditions_From_Options.F90 index ae6f1fc851..abac22bd76 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,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,