sturm_liouville.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2023-2026 N. Tancogne-Dejean, C. Joens
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module sturm_liouville_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use grid_oct_m
  use iso_c_binding
  use lalg_basic_oct_m
  use lapl_operator_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use nl_operator_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use solvers_oct_m
  use space_oct_m
 
  implicit none
 
  private
  public ::                             &
    sturm_liouville_t,                  &
    sturm_liouville_init,               &
    sturm_liouville_end,                &
    sturm_liouville_solve_poisson,      &
    sturm_liouville_solve,              &
    sturm_liouville_solve_from_div
 
 
  type :: sturm_liouville_t
    private
    type(grid_t),    pointer, public :: gr => null()     
    real(real64),    pointer         :: rho(:) => null() 
    real(real64), allocatable :: diag_preconditioner(:)
    real(real64), allocatable :: lapl_rho(:)
    real(real64), allocatable :: diag_lapl(:)
    real(real64), allocatable :: sl_tmp_prod(:)
    real(real64), allocatable :: sl_tmp_lapl(:)
    integer                   :: max_iter = 500
    real(real64)              :: threshold = 1e-7_real64 
    real(real64)              :: inverse_tol = m_epsilon 
    logical                   :: initialized = .false.
  end type sturm_liouville_t
 
contains
 
  ! -------------------------------------------------------------------------------------
  subroutine sturm_liouville_init(this, namespace, gr, space, max_iter, thr, inverse_tol)
    type(sturm_liouville_t),    intent(inout) :: this
    type(namespace_t),          intent(in)    :: namespace
    type(grid_t), target,       intent(in)    :: gr
    type(space_t),              intent(in)    :: space
    integer,          optional, intent(in)    :: max_iter
    real(real64),     optional, intent(in)    :: thr
    real(real64),     optional, intent(in)    :: inverse_tol
 
    character(len=80) :: name
    type(nl_operator_t) :: op(1)
 
    push_sub(sturm_liouville_init)
    assert(.not. this%initialized)
 
    this%gr => gr
    if (present(max_iter)) this%max_iter = max_iter
    if (present(thr))      this%threshold = thr
    if (gr%box%dim == 1)     this%inverse_tol = 1e-20_real64
    if (present(inverse_tol)) this%inverse_tol = inverse_tol
 
    ! Build the Laplacian stencil to extract the diagonal for the Jacobi preconditioner
    name = 'SL preconditioner'
    call derivatives_get_lapl(gr%der, namespace, op, space, name, 1)
    safe_allocate(this%diag_lapl(1:op(1)%np))
    call dnl_operator_operate_diag(op(1), this%diag_lapl)
    call nl_operator_end(op(1))
 
    safe_allocate(this%lapl_rho(1:this%gr%np))
    safe_allocate(this%sl_tmp_prod(1:this%gr%np_part))
    safe_allocate(this%sl_tmp_lapl(1:this%gr%np))
    safe_allocate(this%diag_preconditioner(1:this%gr%np))
 
    this%initialized = .true.
 
    pop_sub(sturm_liouville_init)
  end subroutine sturm_liouville_init
 
 
  ! -------------------------------------------------------------------------------------
  subroutine sturm_liouville_end(this)
    type(sturm_liouville_t), intent(inout) :: this
 
    push_sub(sturm_liouville_end)
 
    safe_deallocate_a(this%diag_lapl)
    safe_deallocate_a(this%diag_preconditioner)
    safe_deallocate_a(this%lapl_rho)
    safe_deallocate_a(this%sl_tmp_prod)
    safe_deallocate_a(this%sl_tmp_lapl)
    nullify(this%gr)
    nullify(this%rho)
    this%initialized = .false.
 
    pop_sub(sturm_liouville_end)
  end subroutine sturm_liouville_end
 
  ! -------------------------------------------------------------------------------------
  subroutine sturm_liouville_solve_poisson(this, psolver, namespace, f, v, with_inverse_density, rho)
    type(sturm_liouville_t), target,  intent(inout) :: this
    type(poisson_t),                  intent(in)    :: psolver
    type(namespace_t),                intent(in)    :: namespace
    real(real64),                     intent(in)    :: f(:,:)
    real(real64),                     intent(inout) :: v(:)
    logical, optional,                intent(in)    :: with_inverse_density
    real(real64), optional, target,   intent(in)    :: rho(:)
 
    real(real64), allocatable :: poisson_rhs(:), tmp_f(:, :)
    real(real64) :: res
    integer :: iter, ip, idm
 
    push_sub(sturm_liouville_solve_poisson)
    assert(this%initialized)
 
    v = m_zero
    if (present(rho)) this%rho => rho
 
    safe_allocate(tmp_f(1:this%gr%np_part, 1:this%gr%box%dim))
    ! We allow the force density F and \f$\frac{F}{\rho}\f$ as initial values of f(:, :)
    if (optional_default(with_inverse_density, .true.)) then
      assert(associated(this%rho))
      do ip = 1, this%gr%np
        do idm = 1, this%gr%box%dim
          tmp_f(ip, idm) = f(ip, idm) / (safe_tol(this%rho(ip), this%inverse_tol))
        end do
      end do
    else
      call lalg_copy(this%gr%np, this%gr%box%dim, f, tmp_f)
    end if
 
    ! \f$ poisson_rhs = \nabla \frac{\vec(F)}{\rho} \f$
    safe_allocate(poisson_rhs(1:this%gr%np))
    call dderivatives_div(this%gr%der, tmp_f(:, :), poisson_rhs(:))
    safe_deallocate_a(tmp_f)
 
    !! Solve the Poisson equation for the initial guess
    !! \f$ v = -\int frac{ \nabla'\cdot\frac{\vec{F} }{\rho}(\vec{r'}){\abs{r-r'}} d^\text{dim}\vec{r'}\f$
    !! For 3D we use the default Poisson solver, for 1D and 2D we use the conjugate gradient solver.
    select case(this%gr%box%dim)
    case(3)
      call dpoisson_solve(psolver, namespace, v, poisson_rhs, all_nodes=.true.)
      call lalg_scal(this%gr%np, -m_fourth / m_pi, v)
    case default
      ! For 1D and 2D we solve the Poisson equation using the conjugate gradient method
      ! Importantly, we cannot use the default Poisson solver,
      ! as we do not want a soft Coulomb here but the Poisson eqaution in 1D.
      iter = 1500
      ! Initialize the mesh_aux pointer needed by dmf_dotu_aux / dmf_nrm2_aux
      call mesh_init_mesh_aux(this%gr)
      call dconjugate_gradients(this%gr%np, v, poisson_rhs, dlaplacian_op, dmf_dotp_aux, iter, res, &
        this%threshold, userdata=[c_loc(this%gr%der)])
      write(message(1),'(a,i6,a)')  "Info: CG converged with ", iter, " iterations."
      write(message(2),'(a,e14.6)')    "Info: The residue is ", res
      call messages_info(2, namespace=namespace)
      ! The dlaplacian_op applies a negative Laplacian.
      ! Hence, we need a negative sign here.
      call lalg_scal(this%gr%np, -m_one, v)
    end select
 
    safe_deallocate_a(poisson_rhs)
    pop_sub(sturm_liouville_solve_poisson)
  end subroutine sturm_liouville_solve_poisson
 
  ! -------------------------------------------------------------------------------------
  subroutine sturm_liouville_solve(this, namespace, rho, rhs, v)
    type(sturm_liouville_t), target,  intent(inout) :: this
    type(namespace_t),                intent(in)    :: namespace
    real(real64), contiguous, target, intent(in)    :: rho(:)
    real(real64), contiguous,         intent(in)    :: rhs(:)
    real(real64), contiguous,         intent(inout) :: v(:)
 
    integer :: iter, ip
    real(real64) :: res
 
    push_sub(sturm_liouville_solve)
 
    assert(this%initialized)
    assert(ubound(rho, dim=1) >= this%gr%np_part)
    assert(ubound(rhs, dim=1) >= this%gr%np)
    assert(ubound(v, dim=1)   >= this%gr%np)
 
    ! Set the density pointer, precompute its Laplacian and set diag_preconditioner for the Jacobi preconditioner
    this%rho => rho
    call mesh_init_mesh_aux(this%gr)
    call dderivatives_lapl(this%gr%der, this%rho, this%lapl_rho)
    !$omp parallel do
    do ip = 1, this%gr%np
      this%diag_preconditioner(ip) = 1 / (safe_tol(this%rho(ip)* this%diag_lapl(ip), this%inverse_tol))
    end do
 
    iter = this%max_iter
    call dqmr_sym_gen_dotu(this%gr%np, v, rhs, sl_operator, dmf_dotu_aux, &
      dmf_nrm2_aux, jacobi_preconditioner, iter, userdata=[c_loc(this)], residue = res, &
      threshold = this%threshold, showprogress = .false.)
 
    write(message(1), '(a, i6, a)') "Info: Sturm-Liouville solver converged in  ", iter, " iterations."
    write(message(2), '(a, es14.6)') "Info: The residue is ", res
    call messages_info(2, namespace=namespace)
 
    ! Cleanup per-solve data
    nullify(this%rho)
 
    pop_sub(sturm_liouville_solve)
  end subroutine sturm_liouville_solve
 
 
  ! -------------------------------------------------------------------------------------
  subroutine sturm_liouville_solve_from_div(this, namespace, rho, f, v, psolver)
    type(sturm_liouville_t), target,  intent(inout) :: this
    type(namespace_t),                intent(in)    :: namespace
    real(real64), target, contiguous, intent(in)    :: rho(:)
    real(real64), contiguous,         intent(in)    :: f(:,:)
    real(real64), contiguous,         intent(inout) :: v(:)
    type(poisson_t), optional,        intent(in)    :: psolver
 
    real(real64), allocatable :: rhs(:), neg_f(:,:)
 
    push_sub(sturm_liouville_solve_from_div)
    assert(this%initialized)
    assert(ubound(rho, dim=1) >= this%gr%np_part)
    assert(ubound(f, dim=1) >= this%gr%np_part)
    assert(ubound(f, dim=2) >= this%gr%box%dim)
    assert(ubound(v, dim=1) >= this%gr%np)
 
    safe_allocate(neg_f(1:this%gr%np_part, 1:this%gr%box%dim))
    neg_f = -f
    ! Initialize the guess with the Poisson solution
    if (present(psolver)) then
      call sturm_liouville_solve_poisson(this, psolver, namespace, neg_f, v, rho=rho)
    end if
 
    safe_allocate(rhs(1:this%gr%np))
    call dderivatives_div(this%gr%der, neg_f, rhs)
    call sturm_liouville_solve(this, namespace, rho, rhs, v)
 
    safe_deallocate_a(rhs)
    safe_deallocate_a(neg_f)
 
    pop_sub(sturm_liouville_solve_from_div)
  end subroutine sturm_liouville_solve_from_div
 
 
  ! -------------------------------------------------------------------------------------
  subroutine sl_operator(x, hx, userdata)
    real(real64), contiguous, intent(in)  :: x(:)
    real(real64), contiguous, intent(out) :: hx(:)
    type(c_ptr),              intent(in)  :: userdata(:)
 
    integer :: ip
    real(real64), allocatable :: v(:), prod(:), lapl_v(:), lapl_product(:)
    type(sturm_liouville_t), pointer :: sl
 
    assert(size(userdata) == 1)
    assert(c_associated(userdata(1)))
    call c_f_pointer(userdata(1), sl)
 
    safe_allocate(v(1:sl%gr%np_part))
    safe_allocate(lapl_v(1:sl%gr%np))
    call lalg_copy(sl%gr%np, x, v)
    call dderivatives_lapl(sl%gr%der, v, lapl_v)
 
    safe_allocate(prod(1:sl%gr%np_part))
    safe_allocate(lapl_product(1:sl%gr%np))
    do ip = 1, sl%gr%np_part
      prod(ip) = sl%rho(ip) * v(ip)
    end do
    call dderivatives_lapl(sl%gr%der, prod, lapl_product, set_bc=.false.)
 
    do ip = 1, sl%gr%np
      hx(ip) = m_half * (sl%rho(ip) * lapl_v(ip) - v(ip) * sl%lapl_rho(ip) + lapl_product(ip))
    end do
 
    safe_deallocate_a(v)
    safe_deallocate_a(prod)
    safe_deallocate_a(lapl_v)
    safe_deallocate_a(lapl_product)
 
  end subroutine sl_operator
 
  subroutine jacobi_preconditioner(x, hx, userdata)
    real(real64), contiguous, intent(in)  :: x(:)
    real(real64), contiguous, intent(out) :: hx(:)
    type(c_ptr),              intent(in)  :: userdata(:)
 
    integer :: ip
    type(sturm_liouville_t), pointer :: sl
 
    assert(size(userdata) == 1)
    assert(c_associated(userdata(1)))
    call c_f_pointer(userdata(1), sl)
 
    !$omp parallel do
    do ip = 1, sl%gr%np
      hx(ip) = x(ip) * sl%diag_preconditioner(ip)
    end do
 
  end subroutine jacobi_preconditioner
 
end module sturm_liouville_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: