scalar_pnpn.f90

Go to the documentation of this file.

! Copyright (c) 2022-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 scalar_pnpn
  use num_types, only : rp, dp
  use rhs_maker, only : rhs_maker_bdf_t, rhs_maker_ext_t, rhs_maker_oifs_t, &
       rhs_maker_ext_fctry, rhs_maker_bdf_fctry, rhs_maker_oifs_fctry
  use scalar_scheme, only : scalar_scheme_t
  use checkpoint_payload, only : checkpoint_payload_t
  use checkpoint, only : chkp_t
  use field, only : field_t
  use scalar_bc_projector, only : scalar_bc_projector_t
  use mesh, only : mesh_t
  use coefs, only : coef_t
  use device, only : host_to_device, device_memcpy, glb_cmd_event, &
       device_event_sync
  use gather_scatter, only : gs_t, gs_op_add, gs_op_min, gs_op_max
  use scalar_residual, only : scalar_residual_t, scalar_residual_factory
  use ax_product, only : ax_t, ax_helm_allocator
  use ax_helm_svv, only : ax_helm_svv_t
  use field_series, only : field_series_t
  use facet_normal, only : facet_normal_t
  use krylov, only : ksp_monitor_t
  use device_math, only : device_add2s2, device_col2
  use time_scheme_controller, only : time_scheme_controller_t
  use projection, only : projection_t
  use math, only : glsc2, col2, add2s2
  use field_math, only : field_col2, field_col3
  use scratch_registry, only : neko_scratch_registry
  use logger, only : neko_log, log_size, neko_log_debug
  use advection, only : advection_t, advection_factory
  use profiler, only : profiler_start_region, profiler_end_region
  use json_utils, only : json_get, json_get_or_default, json_extract_item
  use json_module, only : json_file, json_core, json_value
  use user_intf, only : user_t
  use neko_config, only : neko_bcknd_device
  use time_step_controller, only : time_step_controller_t
  use time_state, only : time_state_t
  use utils, only : neko_error
  use bc, only : bc_t, bc_dirichlet
  use utils, only : neko_varname_len
  use comm, only : neko_comm
  use mpi_f08, only : mpi_allreduce, mpi_integer, mpi_max
  implicit none
  private
 
 
  type, public, extends(scalar_scheme_t) :: scalar_pnpn_t
 
     type(field_t) :: s_res
 
     type(field_t) :: ds
 
     class(ax_t), allocatable :: ax
 
     type(projection_t) :: proj_s
 
     type(scalar_bc_projector_t) :: bc_projector
 
     class(advection_t), allocatable :: adv
 
     ! Time interpolation scheme
     logical :: oifs
 
     ! Advection terms for the oifs method
     type(field_t) :: advs
 
     class(scalar_residual_t), allocatable :: res
 
     class(rhs_maker_ext_t), allocatable :: makeext
 
     class(rhs_maker_bdf_t), allocatable :: makebdf
 
     class(rhs_maker_oifs_t), allocatable :: makeoifs
 
     type(field_t) :: abx1, abx2
 
     type(field_series_t), pointer :: ulag => null()
     type(field_series_t), pointer :: vlag => null()
     type(field_series_t), pointer :: wlag => null()
 
   contains
     procedure, pass(this) :: init => scalar_pnpn_init
     procedure, pass(this) :: restart => scalar_pnpn_restart
     procedure, pass(this) :: register_checkpoint => &
          scalar_pnpn_register_checkpoint
     procedure, pass(this) :: free => scalar_pnpn_free
     procedure, pass(this) :: step => scalar_pnpn_step
     procedure, pass(this) :: apply_strong_bcs => scalar_scheme_apply_strong_bcs
     procedure, pass(this) :: setup_bcs_ => scalar_pnpn_setup_bcs_
  end type scalar_pnpn_t
 
  interface
 
     module subroutine scalar_pnpn_bc_factory(object, scheme, json, coef, user)
       class(bc_t), pointer, intent(inout) :: object
       type(scalar_pnpn_t), intent(in) :: scheme
       type(json_file), intent(inout) :: json
       type(coef_t), target, intent(in) :: coef
       type(user_t), intent(in) :: user
     end subroutine scalar_pnpn_bc_factory
  end interface
 
  interface
 
     module subroutine scalar_pnpn_bc_allocator(object, type_name)
       class(bc_t), pointer, intent(inout) :: object
       character(len=*), intent(in) :: type_name
     end subroutine scalar_pnpn_bc_allocator
  end interface
 
  !
  ! Machinery for injecting user-defined types
  !
 
  abstract interface
     subroutine scalar_pnpn_bc_allocate(obj)
       import bc_t
       class(bc_t), pointer, intent(inout) :: obj
     end subroutine scalar_pnpn_bc_allocate
  end interface
 
  interface
 
     module subroutine register_scalar_pnpn_bc(type_name, allocator)
       character(len=*), intent(in) :: type_name
       procedure(scalar_pnpn_bc_allocate), pointer, intent(in) :: allocator
     end subroutine register_scalar_pnpn_bc
  end interface
 
  type scalar_pnpn_bc_allocator_entry
     character(len=NEKO_VARNAME_LEN) :: type_name
     procedure(scalar_pnpn_bc_allocate), pointer, nopass :: allocator => null()
  end type scalar_pnpn_bc_allocator_entry
 
  type(scalar_pnpn_bc_allocator_entry), allocatable, private :: &
       scalar_pnpn_bc_registry(:)
 
  integer, private :: scalar_pnpn_bc_registry_size = 0
 
  public :: scalar_pnpn_bc_allocator, scalar_pnpn_bc_allocate, &
       register_scalar_pnpn_bc
 
contains
 
  subroutine scalar_pnpn_init(this, msh, coef, gs, params, numerics_params, &
       user, chkp, ulag, vlag, wlag, time_scheme, rho)
    class(scalar_pnpn_t), target, intent(inout) :: this
    type(mesh_t), target, intent(in) :: msh
    type(coef_t), target, intent(in) :: coef
    type(gs_t), target, intent(inout) :: gs
    type(json_file), target, intent(inout) :: params
    type(json_file), target, intent(inout) :: numerics_params
    type(user_t), target, intent(in) :: user
    type(chkp_t), target, intent(inout) :: chkp
    type(field_series_t), target, intent(in) :: ulag, vlag, wlag
    type(time_scheme_controller_t), target, intent(in) :: time_scheme
    type(field_t), target, intent(in) :: rho
    integer :: i
    class(bc_t), pointer :: bc_i
    character(len=15), parameter :: scheme = 'Modular (Pn/Pn)'
    logical :: advection
    real(kind=dp), pointer :: tlag(:), dtlag(:)
 
    call this%free()
 
    ! Initiliaze base type.
    call this%scheme_init(msh, coef, gs, params, scheme, user, rho)
 
    ! Setup backend dependent Ax routines
    if (this%svv_enabled) then
       call ax_helm_allocator(this%ax, type_name = "standard_svv")
       select type (operator => this%ax)
       class is (ax_helm_svv_t)
          operator%svv => this%svv
       end select
    else
       call ax_helm_allocator(this%ax, type_name = "standard")
    end if
 
    ! Setup backend dependent scalar residual routines
    call scalar_residual_factory(this%res)
 
    ! Setup backend dependent summation of extrapolation scheme
    call rhs_maker_ext_fctry(this%makeext)
 
    ! Setup backend dependent contributions to F from lagged BD terms
    call rhs_maker_bdf_fctry(this%makebdf)
 
    ! Setup backend dependent contributions of the OIFS scheme
    call rhs_maker_oifs_fctry(this%makeoifs)
 
    ! Initialize variables specific to this plan
    associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
         dm_xh => this%dm_Xh, nelv => this%msh%nelv)
 
      call this%s_res%init(dm_xh, "s_res")
 
      call this%abx1%init(dm_xh, trim(this%name) // "_abx1")
 
      call this%abx2%init(dm_xh, trim(this%name) // "_abx2")
 
      call this%advs%init(dm_xh, "advs")
 
      call this%ds%init(dm_xh, 'ds')
 
    end associate
 
    ! Set up boundary conditions
    call this%setup_bcs_(user)
 
    do i = 1, this%bcs%size()
       if (this%bcs%bc_type(i) .eq. bc_dirichlet) then
          bc_i => this%bcs%get(i)
          call this%bc_projector%mark(bc_i)
       end if
    end do
 
    ! Initialize projection space
    call this%proj_s%init(this%dm_Xh%size(), this%projection_dim, &
         this%projection_activ_step)
 
    ! Determine the time-interpolation scheme
    call json_get_or_default(numerics_params, 'oifs', this%oifs, .false.)
    ! Point to case checkpoint
    this%chkp => chkp
    ! Initialize advection factory
    call json_get_or_default(params, 'advection', advection, .true.)
    ! OIFS integrates the advection term. With advection disabled, fall back to
    ! the standard BDF history assembly.
    this%oifs = this%oifs .and. advection
 
    this%ulag => ulag
    this%vlag => vlag
    this%wlag => wlag
 
    call chkp%get_time_history(tlag, dtlag)
    call advection_factory(this%adv, numerics_params, this%c_Xh, &
         ulag, vlag, wlag, dtlag, &
         tlag, time_scheme, .not. advection, &
         this%slag)
  end subroutine scalar_pnpn_init
 
  subroutine scalar_pnpn_register_checkpoint(this, chkp)
    class(scalar_pnpn_t), target, intent(inout) :: this
    type(chkp_t), intent(inout) :: chkp
    type(checkpoint_payload_t), pointer :: payload
 
    payload => chkp%add_payload("scalars/" // trim(this%name))
    call payload%add_field(this%s)
    call payload%add_series(this%slag)
    call payload%add_field(this%abx1)
    call payload%add_field(this%abx2)
 
  end subroutine scalar_pnpn_register_checkpoint
 
  ! Restarts the scalar from a checkpoint
  subroutine scalar_pnpn_restart(this, chkp)
    class(scalar_pnpn_t), target, intent(inout) :: this
    type(chkp_t), intent(inout) :: chkp
    integer :: i, n
    class(bc_t), pointer :: bc_i
    logical :: interpolated
 
    n = this%s%dof%size()
 
    ! Lag fields are restored through the checkpoint's fsp mechanism
 
    ! The restored fields are continuous unless the checkpoint was written
    ! on another mesh or at another polynomial order and was interpolated
    ! on the way in. Only then do the copies of a node shared between
    ! elements need averaging: scale by the inverse multiplicity, then sum
    ! the copies with a gather-scatter. On a plain restart that is the
    ! identity in exact arithmetic but not in floating point, and would put
    ! about one ulp of error on every shared node. Same guard as the fluid.
    interpolated = allocated(chkp%previous_mesh%elements) .or. &
         chkp%previous_Xh%lx .ne. this%Xh%lx
 
    if (interpolated) then
       call col2(this%s%x, this%c_Xh%mult, n)
       call col2(this%slag%lf(1)%x, this%c_Xh%mult, n)
       call col2(this%slag%lf(2)%x, this%c_Xh%mult, n)
    end if
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%s%x, this%s%x_d, &
            n, host_to_device, sync = .false.)
       call device_memcpy(this%slag%lf(1)%x, this%slag%lf(1)%x_d, &
            n, host_to_device, sync = .false.)
       call device_memcpy(this%slag%lf(2)%x, this%slag%lf(2)%x_d, &
            n, host_to_device, sync = .false.)
       call device_memcpy(this%abx1%x, this%abx1%x_d, &
            n, host_to_device, sync = .false.)
       call device_memcpy(this%abx2%x, this%abx2%x_d, &
            n, host_to_device, sync = .false.)
       call device_memcpy(this%advs%x, this%advs%x_d, &
            n, host_to_device, sync = .false.)
    end if
 
    if (interpolated) then
       call this%gs_Xh%op(this%s, gs_op_add)
       call this%gs_Xh%op(this%slag%lf(1), gs_op_add)
       call this%gs_Xh%op(this%slag%lf(2), gs_op_add)
    end if
 
    ! Restore scalar bcs that need it. This is a no op in most bcs.
    do i = 1, this%bcs%size()
       bc_i => this%bcs%get(i)
       call bc_i%restart(this%s, this%slag)
    end do
 
    nullify(bc_i)
 
  end subroutine scalar_pnpn_restart
 
  subroutine scalar_pnpn_free(this)
    class(scalar_pnpn_t), intent(inout) :: this
 
    ! Release operator references before scheme_free deallocates their targets.
    if (allocated(this%Ax)) then
       call this%Ax%free()
       deallocate(this%Ax)
    end if
 
    !Deallocate scalar field
    call this%scheme_free()
 
    call this%bc_projector%free()
    call this%proj_s%free()
 
    call this%s_res%free()
 
    call this%ds%free()
 
    call this%abx1%free()
    call this%abx2%free()
 
    call this%advs%free()
 
    if (allocated(this%adv)) then
       call this%adv%free()
       deallocate(this%adv)
    end if
 
    nullify(this%ulag)
    nullify(this%vlag)
    nullify(this%wlag)
 
    if (allocated(this%res)) then
       deallocate(this%res)
    end if
 
    if (allocated(this%makeext)) then
       deallocate(this%makeext)
    end if
 
    if (allocated(this%makebdf)) then
       deallocate(this%makebdf)
    end if
 
    if (allocated(this%makeoifs)) then
       deallocate(this%makeoifs)
    end if
 
  end subroutine scalar_pnpn_free
 
  subroutine scalar_pnpn_step(this, time, ext_bdf, dt_controller, &
       ksp_results)
    class(scalar_pnpn_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
    type(time_scheme_controller_t), intent(in) :: ext_bdf
    type(time_step_controller_t), intent(in) :: dt_controller
    type(ksp_monitor_t), intent(inout) :: ksp_results
    type(field_t), pointer :: rho_cp
    integer :: rho_cp_index
    ! Number of degrees of freedom
    integer :: n
 
    if (this%freeze) return
 
    n = this%dm_Xh%size()
    call neko_scratch_registry%request_field(rho_cp, rho_cp_index, .false.)
 
    call profiler_start_region(trim(this%name), 2)
    associate(u => this%u, v => this%v, w => this%w, s => this%s, &
         cp => this%cp, rho => this%rho, lambda_tot => this%lambda_tot, &
         ds => this%ds, &
         s_res => this%s_res, &
         ax => this%Ax, f_xh => this%f_Xh, xh => this%Xh, &
         c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
         slag => this%slag, oifs => this%oifs, &
         projection_dim => this%projection_dim, &
         msh => this%msh, res => this%res, makeoifs => this%makeoifs, &
         makeext => this%makeext, makebdf => this%makebdf, &
         t => time%t, tstep => time%tstep, dt => time%dt)
 
      ! Logs extra information the log level is NEKO_LOG_DEBUG or above.
      call print_debug(this)
 
      ! Update material properties and their pointwise product.
      call this%update_material_properties(time)
      call field_col3(rho_cp, rho, cp, n)
 
      ! Update the SVV coefficient if SVV is enabled.
      if (this%svv_enabled) then
         call this%svv%update(rho_cp, tstep)
      end if
 
      ! Compute the source terms
      call this%source_term%compute(time)
 
      if (oifs) then
         ! The fluid step has already advanced u, v, and w to the new time.
         ! Its first lag fields contain the velocity at tlag(1), which is the
         ! latest time represented by the OIFS interpolation history.
         call this%adv%compute_scalar(this%ulag%lf(1), this%vlag%lf(1), &
              this%wlag%lf(1), s, this%advs, &
              xh, this%c_Xh, dm_xh%size())
      else
         ! Add the advection operators to the right-hand side.
         call this%adv%compute_scalar(u, v, w, s, f_xh, &
              xh, this%c_Xh, dm_xh%size())
      end if
 
      ! Scale the volumetric source and advection terms by rho * cp.
      call field_col2(f_xh, rho_cp, n)
 
      ! Add weak boundary fluxes without scaling them by rho * cp.
      call this%bcs%apply_scalar(f_xh%x, n, time, .false.)
 
      ! Extrapolate the already scaled explicit right-hand side.
      call makeext%compute_scalar(this%abx1, this%abx2, f_xh%x, &
           ext_bdf%advection_coeffs%x, n)
 
      if (oifs) then
         call makeoifs%compute_scalar(this%advs%x, f_xh%x, &
              rho_cp, real(dt, kind=rp), n)
      else
 
         ! Add the RHS contributions coming from the BDF scheme.
         call makebdf%compute_scalar(slag, f_xh%x, s, c_xh%B, &
              rho_cp, real(dt, kind=rp), ext_bdf%diffusion_coeffs%x, &
              ext_bdf%ndiff, n)
      end if
 
      call slag%update()
 
      call this%apply_strong_bcs(time)
 
      ! Compute scalar residual.
      call profiler_start_region(trim(this%name) // '_residual', 20)
      call res%compute(ax, s, s_res, f_xh, c_xh, msh, xh, lambda_tot, &
           rho_cp, ext_bdf%diffusion_coeffs%x(1), &
           real(dt, kind=rp), dm_xh%size())
 
      call gs_xh%op(s_res, gs_op_add)
 
      ! Zero-out residual at Dirichlet nodes before solving.
      call this%bc_projector%apply(s_res%x, dm_xh%size())
 
      call profiler_end_region(trim(this%name) // '_residual', 20)
 
      call this%proj_s%pre_solving(s_res%x, tstep, c_xh, n, dt_controller)
 
      call this%pc%update()
      call profiler_start_region(trim(this%name) // '_solve', 21)
      ksp_results = this%ksp%solve(ax, ds, s_res%x, n, &
           c_xh, this%bc_projector, gs_xh)
      ksp_results%name = trim(this%name)
      call profiler_end_region(trim(this%name) // '_solve', 21)
 
      call this%proj_s%post_solving(ds%x, ax, c_xh, this%bc_projector, gs_xh, &
           n, tstep, dt_controller)
 
      ! Update the solution
      if (neko_bcknd_device .eq. 1) then
         call device_add2s2(s%x_d, ds%x_d, 1.0_rp, n)
      else
         call add2s2(s%x, ds%x, 1.0_rp, n)
      end if
 
    end associate
    call neko_scratch_registry%relinquish_field(rho_cp_index)
    call profiler_end_region(trim(this%name), 2)
  end subroutine scalar_pnpn_step
 
  subroutine print_debug(this)
    class(scalar_pnpn_t), intent(inout) :: this
    character(len=LOG_SIZE) :: log_buf
    integer :: n
 
    n = this%dm_Xh%size()
 
    write(log_buf, '(A, A, E15.7, A, E15.7, A, E15.7)') 'Scalar debug', &
         ' l2norm s', glsc2(this%s%x, this%s%x, n), &
         ' slag1', glsc2(this%slag%lf(1)%x, this%slag%lf(1)%x, n), &
         ' slag2', glsc2(this%slag%lf(2)%x, this%slag%lf(2)%x, n)
    call neko_log%message(log_buf, lvl = neko_log_debug)
    write(log_buf, '(A, A, E15.7, A, E15.7)') 'Scalar debug2', &
         ' l2norm abx1', glsc2(this%abx1%x, this%abx1%x, n), &
         ' abx2', glsc2(this%abx2%x, this%abx2%x, n)
    call neko_log%message(log_buf, lvl = neko_log_debug)
  end subroutine print_debug
 
  subroutine scalar_pnpn_setup_bcs_(this, user)
    class(scalar_pnpn_t), target, intent(inout) :: this
    type(user_t), target, intent(in) :: user
    integer :: i, j, n_bcs, zone_size, global_zone_size, ierr
    type(json_core) :: core
    type(json_value), pointer :: bc_object
    type(json_file) :: bc_subdict
    class(bc_t), pointer :: bc_i
    logical :: found
    ! Monitor which boundary zones have been marked
    logical, allocatable :: marked_zones(:)
    integer, allocatable :: zone_indices(:)
    character(len=256) :: error_msg
 
    if (this%params%valid_path('boundary_conditions')) then
       call this%params%info('boundary_conditions', &
            n_children = n_bcs)
       call this%params%get_core(core)
       call this%params%get('boundary_conditions', bc_object, found)
 
       call this%bcs%init(n_bcs)
 
       allocate(marked_zones(size(this%msh%labeled_zones)))
       marked_zones = .false.
 
       do i = 1, n_bcs
          ! Create a new json containing just the subdict for this bc
          call json_extract_item(core, bc_object, i, bc_subdict)
 
          ! Check that we are not trying to assing a bc to zone, for which one
          ! has already been assigned and that the zone has more than 0 size
          ! in the mesh.
          call json_get(bc_subdict, "zone_indices", zone_indices)
 
          do j = 1, size(zone_indices)
             zone_size = this%msh%labeled_zones(zone_indices(j))%size
             call mpi_allreduce(zone_size, global_zone_size, 1, &
                  mpi_integer, mpi_max, neko_comm, ierr)
 
             if (global_zone_size .eq. 0) then
                write(error_msg, '(A, I0, A, A, I0, A)') &
                     "Zone index ", zone_indices(j), &
                     " is invalid as this zone has 0 size, meaning it ", &
                     "does not exist in the mesh. Check scalar boundary ", &
                     "condition ", i, "."
                call neko_error(error_msg)
                error stop
             end if
 
             if (marked_zones(zone_indices(j))) then
                write(error_msg, '(A, I0, A, A, A, A)')&
                     "Zone with index ", zone_indices(j), &
                     " has already been assigned a boundary condition. ", &
                     "Please check your boundary_conditions entry for the ", &
                     "scalar and make sure that each zone index appears only ",&
                     "in a single boundary condition."
                call neko_error(error_msg)
                error stop
             else
                marked_zones(zone_indices(j)) = .true.
             end if
          end do
 
          bc_i => null()
 
          call scalar_pnpn_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
          call this%bcs%append(bc_i)
       end do
 
       ! Make sure all labeled zones with non-zero size have been marked
       do i = 1, size(this%msh%labeled_zones)
          if ((this%msh%labeled_zones(i)%size .gt. 0) .and. &
               (.not. marked_zones(i))) then
             write(error_msg, '(A, I0)') &
                  "No scalar boundary condition assigned to zone ", i
             call neko_error(error_msg)
             error stop
          end if
       end do
    else
       ! Check that there are no labeled zones, i.e. all are periodic.
       do i = 1, size(this%msh%labeled_zones)
          if (this%msh%labeled_zones(i)%size .gt. 0) then
             write(error_msg, '(A, A)') &
                  "No boundary_conditions entry in the case file for scalar ", &
                  this%s%name
             call neko_error(error_msg)
             error stop
          end if
       end do
 
       ! For a pure periodic case, we still need to initilise the bc lists
       ! to a zero size to avoid issues with apply() in step()
       call this%bcs%init()
 
    end if
  end subroutine scalar_pnpn_setup_bcs_
 
  subroutine scalar_scheme_apply_strong_bcs(this, time)
    class(scalar_pnpn_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
 
    integer :: i
    class(bc_t), pointer :: bc_i
    bc_i => null()
 
    ! First apply call, sets the Dirichlet value, let's call it d.
    call this%bcs%apply(this%s, time = time, strong = .true.)
    ! If we now have local nodes sharing the same global node, and with
    ! some nodes not masked as Dirichlet, the node which *is* masked
    ! will have the value d, and the the other ones just some value u.
    ! Take a nodewise minimum between the local nodes.
    ! Now, all local nodes store m = min(d, u)
    call this%gs_Xh%op(this%s, gs_op_min, glb_cmd_event)
    call device_event_sync(glb_cmd_event)
 
    ! Second apply call, so Dirichlet nodes again store d, the rest still store
    ! m, where m < d by construction.
    call this%bcs%apply(this%s, time = time, strong = .true.)
    ! Now apply a max, which guarantees that d wins and gets stored in all the
    ! local nodes.
    call this%gs_Xh%op(this%s, gs_op_max, glb_cmd_event)
    call device_event_sync(glb_cmd_event)
 
    ! Reset updated flags
    do i = 1, this%bcs%size()
       bc_i => this%bcs%get(i)
       bc_i%updated = .false.
    end do
    nullify(bc_i)
 
  end subroutine scalar_scheme_apply_strong_bcs
 
 
end module scalar_pnpn