PDE_filter.f90

Go to the documentation of this file.

! Copyright (c) 2023-2026, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
 
module pde_filter
  use num_types, only : rp
  use json_module, only : json_file
  use registry, only : neko_registry
  use field, only : field_t
  use coefs, only : coef_t
  use ax_product, only : ax_t, ax_helm_allocator
  use krylov, only : ksp_t, ksp_monitor_t, krylov_solver_factory
  use precon, only : pc_t, precon_allocator, precon_destroy
  use bc_list, only : bc_list_t
  use scalar_bc_projector, only : scalar_bc_projector_t
  use neumann, only : neumann_t
  use profiler, only : profiler_start_region, profiler_end_region
  use gather_scatter, only : gs_t, gs_op_add
  use pnpn_residual, only : pnpn_prs_res_t
  use mesh, only : mesh_t, neko_msh_max_zlbls, neko_msh_max_zlbl_len
  use registry, only : neko_registry
  use filter, only : filter_t
  use scratch_registry, only : neko_scratch_registry
  use field_math, only : field_copy, field_add3
  use coefs, only : coef_t
  use logger, only : neko_log, log_size
  use neko_config, only : neko_bcknd_device
  use dofmap, only : dofmap_t
  use jacobi, only : jacobi_t
  use device_jacobi, only : device_jacobi_t
  use sx_jacobi, only : sx_jacobi_t
  use utils, only : neko_error
  use device_math, only : device_cfill, device_subcol3, device_cmult
  use json_utils, only : json_get, json_get_or_default
  implicit none
  private
 
  type, public, extends(filter_t) :: pde_filter_t
 
     class(ax_t), allocatable :: ax
     type(ksp_monitor_t) :: ksp_results(1)
     class(ksp_t), allocatable :: ksp_filt
     class(pc_t), allocatable :: pc_filt
     type(scalar_bc_projector_t) :: bc_projector_filt
 
     ! Inputs from the user
     real(kind=rp) :: r
     real(kind=rp) :: abstol_filt
     integer :: ksp_max_iter
     character(len=:), allocatable :: ksp_solver
     ! > preconditioner type
     character(len=:), allocatable :: precon_type_filt
     integer :: ksp_n, n, i
 
 
 
   contains
     procedure, pass(this) :: init => pde_filter_init_from_json
     procedure, pass(this) :: init_from_components => &
          pde_filter_init_from_components
     procedure, pass(this) :: free => pde_filter_free
     procedure, pass(this) :: apply => pde_filter_apply
  end type pde_filter_t
 
contains
 
  subroutine pde_filter_init_from_json(this, json, coef)
    class(pde_filter_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    type(coef_t), target, intent(in) :: coef
    real(kind=rp) :: r, tol
    integer :: max_iter
    character(len=:), allocatable :: ksp_solver, precon_type
 
    ! user parameters
    call json_get(json, "radius", r)
    call json_get_or_default(json, "tolerance", tol, 1e-10_rp)
    call json_get_or_default(json, "max_iter", max_iter, 200)
    call json_get_or_default(json, "solver", ksp_solver, "cg")
    call json_get_or_default(json, "preconditioner", precon_type, "jacobi")
 
    call this%init_from_components(coef, r, tol, max_iter, ksp_solver, &
         precon_type)
 
  end subroutine pde_filter_init_from_json
 
  subroutine pde_filter_init_from_components(this, coef, r, tol, max_iter, &
       ksp_solver, precon_type)
    class(pde_filter_t), intent(inout) :: this
    type(coef_t), target, intent(in) :: coef
    real(kind=rp), intent(in) :: r, tol
    integer, intent(in) :: max_iter
    character(len=*), intent(in) :: ksp_solver, precon_type
    integer :: n
 
    call this%init_base(coef)
 
    this%r = r
    this%abstol_filt = tol
    this%ksp_max_iter = max_iter
    this%ksp_solver = ksp_solver
    this%precon_type_filt = precon_type
 
    ! set the number of dofs
    n = this%coef%dof%size()
 
    ! Setup backend dependent Ax routines
    call ax_helm_allocator(this%Ax, type_name = "standard")
 
    ! set up krylov solver
    call krylov_solver_factory(this%ksp_filt, n, this%ksp_solver, &
         this%ksp_max_iter, this%abstol_filt)
 
    ! set up preconditioner
    call filter_precon_factory(this%pc_filt, this%ksp_filt, &
         this%coef, this%coef%dof, &
         this%coef%gs_h, &
         this%bc_projector_filt, this%precon_type_filt)
 
  end subroutine pde_filter_init_from_components
 
  subroutine pde_filter_free(this)
    class(pde_filter_t), intent(inout) :: this
 
    if (allocated(this%Ax)) then
       call this%Ax%free()
       deallocate(this%Ax)
    end if
 
    if (allocated(this%ksp_filt)) then
       call this%ksp_filt%free()
       deallocate(this%ksp_filt)
    end if
 
    if (allocated(this%pc_filt)) then
       call precon_destroy(this%pc_filt)
       deallocate(this%pc_filt)
    end if
 
    if (allocated(this%ksp_solver)) then
       deallocate(this%ksp_solver)
    end if
 
    if (allocated(this%precon_type_filt)) then
       deallocate(this%precon_type_filt)
    end if
 
    call this%bc_projector_filt%free()
 
    call this%free_base()
 
  end subroutine pde_filter_free
 
  subroutine pde_filter_apply(this, F_out, F_in)
    class(pde_filter_t), intent(inout) :: this
    type(field_t), intent(in) :: F_in
    type(field_t), intent(inout) :: F_out
    integer :: n, i
    type(field_t), pointer :: RHS, d_F_out
    character(len=LOG_SIZE) :: log_buf
    integer :: temp_indices(2)
 
    n = this%coef%dof%size()
    call neko_scratch_registry%request_field(rhs, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(d_f_out, temp_indices(2), .false.)
    ! in a similar fasion to pressure/velocity, we will solve for d_F_out.
 
    ! to improve convergence, we use F_in as an initial guess for F_out.
    ! so F_out = F_in + d_F_in.
 
    ! Defining the operator A = -r^2 \nabla^2 + I
    ! the system changes from:
    ! A (F_out) = F_in
    ! to
    ! A (d_F_out) = F_in - A(F_in)
 
    ! set up Helmholtz operators and RHS
    if (neko_bcknd_device .eq. 1) then
       call device_cfill(this%coef%h1_d, &
            (this%r / (2.0_rp * sqrt(3.0_rp)))**2, n)
       call device_cfill(this%coef%h2_d, 1.0_rp, n)
    else
       ! h1 is already negative in its definition
       this%coef%h1 = (this%r / (2.0_rp * sqrt(3.0_rp)))**2
       ! ax_helm includes the mass matrix in h2
       this%coef%h2 = 1.0_rp
    end if
    this%coef%ifh2 = .true.
 
    ! compute the A(F_in) component of the RHS
    ! (note, to be safe with the inout intent we first copy F_in to the
    !  temporary d_F_out)
    call field_copy(d_f_out, f_in)
    call this%Ax%compute(rhs%x, d_f_out%x, this%coef, this%coef%msh, &
         this%coef%Xh)
 
    if (neko_bcknd_device .eq. 1) then
       call device_subcol3(rhs%x_d, f_in%x_d, this%coef%B_d, n)
       call device_cmult(rhs%x_d, -1.0_rp, n)
    else
       do i = 1, n
          ! mass matrix should be included here
          rhs%x(i,1,1,1) = f_in%x(i,1,1,1) * this%coef%B(i,1,1,1) &
               - rhs%x(i,1,1,1)
       end do
    end if
 
    ! gather scatter
    call this%coef%gs_h%op(rhs, gs_op_add)
 
    ! set BCs
    call this%bc_projector_filt%apply(rhs%x, n)
 
    ! Solve Helmholtz equation
    call profiler_start_region("filter solve")
    this%ksp_results(1) = &
         this%ksp_filt%solve(this%Ax, d_f_out, rhs%x, n, this%coef, &
         this%bc_projector_filt, this%coef%gs_h)
 
    call profiler_end_region
 
    ! add result
    call field_add3(f_out, f_in, d_f_out)
    ! update preconditioner (needed?)
    call this%pc_filt%update()
 
    ! write it all out
    call neko_log%section('PDE Filter')
 
    write(log_buf, '(A,A,A)') 'Iterations:   ', 'Start residual:     ', &
         'Final residual:'
    call neko_log%message(log_buf)
    write(log_buf, '(I11,3x, E15.7,5x, E15.7)') this%ksp_results%iter, &
         this%ksp_results%res_start, this%ksp_results%res_final
    call neko_log%message(log_buf)
    call neko_log%end_section()
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine pde_filter_apply
 
  subroutine filter_precon_factory(pc, ksp, coef, dof, gs, bc_projector, &
       pctype)
    class(pc_t), allocatable, target, intent(inout) :: pc
    class(ksp_t), target, intent(inout) :: ksp
    type(coef_t), target, intent(in) :: coef
    type(dofmap_t), target, intent(in) :: dof
    type(gs_t), target, intent(inout) :: gs
    type(scalar_bc_projector_t), target, intent(inout) :: bc_projector
    character(len=*) :: pctype
 
    call precon_allocator(pc, pctype)
 
    select type (pcp => pc)
    type is (jacobi_t)
       call pcp%init(coef, dof, gs)
    type is (sx_jacobi_t)
       call pcp%init(coef, dof, gs)
    type is (device_jacobi_t)
       call pcp%init(coef, dof, gs)
    end select
 
    call ksp%set_pc(pc)
 
  end subroutine filter_precon_factory
 
end module pde_filter