Projection for stress formulation

src/.depends

@@ -204,7 +204,7 @@ fluid/stress_formulation/bcknd/cpu/pnpn_res_stress_cpu.lo : fluid/stress_formula
 fluid/stress_formulation/bcknd/device/pnpn_res_stress_device.lo : fluid/stress_formulation/bcknd/device/pnpn_res_stress_device.F90 math/bcknd/device/device_math.lo math/bcknd/device/device_mathops.lo sem/space.lo config/num_types.lo mesh/mesh.lo field/scratch_registry.lo fluid/pnpn_res.lo bc/facet_normal.lo sem/coef.lo math/ax.lo field/field.lo math/operators.lo common/utils.lo gs/gather_scatter.lo 
 fluid/fluid_scheme.lo : fluid/fluid_scheme.f90 common/gradient_jump_penalty.lo bc/shear_stress.lo bc/wall_model_bc.lo math/field_math.lo common/time_step_controller.lo field/field_series.lo common/utils.lo common/user_intf.lo field/scratch_registry.lo common/json_utils.lo field/field_registry.lo common/log.lo math/operators.lo time_schemes/time_scheme_controller.lo math/bcknd/device/device_math.lo math/math.lo mesh/mesh.lo bc/bc.lo fluid/fluid_stats.lo krylov/precon.lo multigrid/phmg.lo krylov/pc_hsmg.lo krylov/bcknd/device/pc_jacobi_device.lo krylov/bcknd/sx/pc_jacobi_sx.lo krylov/bcknd/cpu/pc_jacobi.lo bc/field_dirichlet_vector.lo bc/field_dirichlet.lo bc/symmetry.lo bc/dong_outflow.lo bc/dirichlet.lo bc/blasius.lo bc/usr_inflow.lo bc/inflow.lo bc/wall.lo sem/coef.lo krylov/krylov.lo sem/dofmap.lo sem/space.lo field/field.lo fluid/fluid_source_term.lo comm/comm.lo config/num_types.lo fluid/mean_flow.lo common/checkpoint.lo config/neko_config.lo fluid/mean_sqr_flow.lo gs/gather_scatter.lo 
 fluid/fluid_aux.lo : fluid/fluid_aux.f90 common/utils.lo krylov/krylov.lo config/num_types.lo common/log.lo 
-fluid/fluid_pnpn.lo : fluid/fluid_pnpn.f90 math/operators.lo math/field_math.lo common/utils.lo bc/bc.lo math/mathops.lo math/math.lo config/neko_config.lo gs/gs_ops.lo common/time_step_controller.lo common/user_intf.lo mesh/mesh.lo comm/comm.lo bc/non_normal.lo bc/facet_normal.lo bc/dirichlet.lo field/field.lo math/ax.lo common/json_utils.lo common/profiler.lo fluid/advection.lo device/device.lo common/projection.lo time_schemes/time_scheme_controller.lo fluid/fluid_aux.lo math/bcknd/device/device_mathops.lo fluid/fluid_scheme.lo fluid/fluid_volflow.lo common/rhs_maker.lo fluid/pnpn_res.lo krylov/krylov.lo config/num_types.lo 
+fluid/fluid_pnpn.lo : fluid/fluid_pnpn.f90 math/operators.lo math/field_math.lo common/utils.lo bc/bc.lo math/mathops.lo math/math.lo config/neko_config.lo gs/gs_ops.lo common/time_step_controller.lo common/user_intf.lo mesh/mesh.lo comm/comm.lo bc/non_normal.lo bc/facet_normal.lo bc/dirichlet.lo field/field.lo math/ax.lo common/json_utils.lo common/profiler.lo fluid/advection.lo device/device.lo common/projection_vel.lo common/projection.lo time_schemes/time_scheme_controller.lo fluid/fluid_aux.lo math/bcknd/device/device_mathops.lo fluid/fluid_scheme.lo fluid/fluid_volflow.lo common/rhs_maker.lo fluid/pnpn_res.lo krylov/krylov.lo config/num_types.lo 
 fluid/fluid_fctry.lo : fluid/fluid_fctry.f90 common/utils.lo fluid/fluid_pnpn.lo fluid/fluid_scheme.lo 
 fluid/fluid_volflow.lo : fluid/fluid_volflow.f90 math/ax.lo bc/bc.lo field/scratch_registry.lo common/json_utils.lo gs/gather_scatter.lo math/bcknd/device/device_mathops.lo math/bcknd/device/device_math.lo config/neko_config.lo comm/comm.lo math/math.lo time_schemes/time_scheme_controller.lo sem/coef.lo field/field.lo sem/dofmap.lo krylov/precon.lo krylov/krylov.lo math/mathops.lo config/num_types.lo math/operators.lo 
 fluid/pnpn_res.lo : fluid/pnpn_res.f90 config/num_types.lo mesh/mesh.lo sem/space.lo bc/facet_normal.lo sem/coef.lo field/field.lo math/ax.lo gs/gather_scatter.lo 
@@ -239,6 +239,7 @@ math/bcknd/xsmm/ax_helm_xsmm.lo : math/bcknd/xsmm/ax_helm_xsmm.F90 math/mxm_wrap
 math/bcknd/device/ax_helm_device.lo : math/bcknd/device/ax_helm_device.F90 device/device.lo math/bcknd/device/device_math.lo mesh/mesh.lo sem/space.lo sem/coef.lo config/num_types.lo math/ax_helm.lo 
 math/bcknd/device/ax_helm_full_device.lo : math/bcknd/device/ax_helm_full_device.F90 common/utils.lo device/device.lo math/bcknd/device/device_math.lo mesh/mesh.lo sem/space.lo sem/coef.lo config/num_types.lo math/ax_helm_full.lo 
 common/projection.lo : common/projection.f90 common/time_step_controller.lo common/utils.lo common/log.lo common/profiler.lo common/bcknd/device/device_projection.lo math/bcknd/device/device_math.lo device/device.lo config/neko_config.lo gs/gather_scatter.lo comm/comm.lo bc/bc.lo math/ax.lo sem/coef.lo math/math.lo config/num_types.lo 
+common/projection_vel.lo : common/projection_vel.f90 common/projection.lo common/time_step_controller.lo common/utils.lo common/log.lo common/profiler.lo common/bcknd/device/device_projection.lo math/bcknd/device/device_math.lo device/device.lo config/neko_config.lo gs/gather_scatter.lo comm/comm.lo bc/bc.lo math/ax.lo sem/coef.lo math/math.lo config/num_types.lo 
 common/bcknd/device/device_projection.lo : common/bcknd/device/device_projection.F90 common/utils.lo config/num_types.lo 
 comm/parmetis.lo : comm/parmetis.F90 mesh/mesh.lo field/mesh_field.lo config/num_types.lo common/utils.lo mesh/point.lo comm/comm.lo 
 comm/redist.lo : comm/redist.f90 mesh/element.lo mesh/facet_zone.lo io/format/nmsh.lo mesh/mesh.lo comm/comm.lo mesh/curve.lo adt/stack.lo mesh/point.lo adt/htable.lo comm/mpi_types.lo field/mesh_field.lo 

src/Makefile.am

@@ -243,6 +243,7 @@ neko_fortran_SOURCES = \
 	math/bcknd/device/ax_helm_device.F90\
 	math/bcknd/device/ax_helm_full_device.F90\
 	common/projection.f90\
+	common/projection_vel.f90\
 	common/bcknd/device/device_projection.F90\
 	comm/parmetis.F90\
 	comm/redist.f90\

src/common/projection.f90

@@ -82,6 +82,7 @@ module projection
 
   implicit none
   private
+  public :: cpu_proj_ortho, device_proj_ortho
 
   type, public ::  projection_t
      real(kind=rp), allocatable :: xx(:,:)

src/common/projection_vel.f90

@@ -0,0 +1,256 @@
+! Copyright (c) 2024, 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.
+!
+!> Project x onto X , the space of old solutions and back again
+!> Couple projections for velocity
+module projection_vel
+  use num_types, only : rp, c_rp
+  use math, only : rzero, glsc3, add2, copy, cmult
+  use coefs, only : coef_t
+  use ax_product, only : ax_t
+  use bc, only : bc_list_t, bc_list_apply_scalar
+  use comm
+  use gather_scatter, only : gs_t, GS_OP_ADD
+  use neko_config, only : NEKO_BCKND_DEVICE
+  use device
+  use device_math, only : device_glsc3, device_add2s2, device_cmult, &
+       device_rzero, device_copy, device_add2, device_add2s2_many, &
+       device_glsc3_many
+  use device_projection, only : device_proj_on, device_project_ortho
+  use profiler, only : profiler_start_region, profiler_end_region
+  use logger, only : LOG_SIZE, neko_log
+  use utils, only : neko_warning
+  use, intrinsic :: iso_c_binding
+  use time_step_controller, only : time_step_controller_t
+  use projection, only : projection_t, cpu_proj_ortho, device_proj_ortho
+
+  implicit none
+  private
+
+  type, public ::  projection_vel_t
+     type(projection_t) :: proj_u, proj_v, proj_w
+     integer :: activ_step
+     integer :: L
+   contains
+     procedure, pass(this) :: init => projection_init_vel
+     procedure, pass(this) :: free => projection_free_vel
+     procedure, pass(this) :: pre_solving => projection_pre_solving_vel
+     procedure, pass(this) :: post_solving => projection_post_solving_vel
+     procedure, pass(this) :: project_back => bcknd_project_back_vel
+  end type projection_vel_t
+
+contains
+
+  subroutine projection_init_vel(this, n, L, activ_step)
+    class(projection_vel_t), target, intent(inout) :: this
+    integer, intent(in) :: n
+    integer, optional, intent(in) :: L, activ_step
+    integer :: i
+    integer(c_size_t) :: ptr_size
+    type(c_ptr) :: ptr
+    real(c_rp) :: dummy
+
+    call this%free()
+
+    call this%proj_u%init(n, L, activ_step)
+    call this%proj_v%init(n, L, activ_step)
+    call this%proj_w%init(n, L, activ_step)
+
+    this%L = L
+    this%activ_step = activ_step
+
+
+  end subroutine projection_init_vel
+
+  subroutine projection_free_vel(this)
+    class(projection_vel_t), intent(inout) :: this
+    integer :: i
+
+    call this%proj_u%free()
+    call this%proj_v%free()
+    call this%proj_w%free()
+
+  end subroutine projection_free_vel
+
+  subroutine projection_pre_solving_vel(this, b_u, b_v, b_w, tstep, coef, n, &
+                                        dt_controller, string)
+    class(projection_vel_t), intent(inout) :: this
+    integer, intent(inout) :: n
+    real(kind=rp), intent(inout), dimension(n) :: b_u, b_v, b_w
+    integer, intent(in) :: tstep
+    class(coef_t), intent(inout) :: coef
+    type(time_step_controller_t), intent(in) :: dt_controller
+    character(len=*), optional :: string
+
+    call this%proj_u%pre_solving(b_u, tstep, coef, n, dt_controller, string)
+    call this%proj_v%pre_solving(b_v, tstep, coef, n, dt_controller, string)
+    call this%proj_w%pre_solving(b_w, tstep, coef, n, dt_controller, string)
+
+  end subroutine projection_pre_solving_vel
+
+  subroutine projection_post_solving_vel(this, x_u, x_v, x_w, Ax, coef, &
+                                         bclst_u, bclst_v, bclst_w, &
+                                         gs_h, n, tstep, dt_controller)
+    class(projection_vel_t), intent(inout) :: this
+    integer, intent(inout) :: n
+    class(Ax_t), intent(inout) :: Ax
+    class(coef_t), intent(inout) :: coef
+    class(bc_list_t), intent(inout) :: bclst_u, bclst_v, bclst_w
+    type(gs_t), intent(inout) :: gs_h
+    real(kind=rp), intent(inout), dimension(n) :: x_u, x_v, x_w
+    integer, intent(in) :: tstep
+    type(time_step_controller_t), intent(in) :: dt_controller
+
+    ! Here we assume the projection space sizes and activate steps
+    ! for all three velocity equations are the same
+    if (tstep .gt. this%activ_step .and. this%L .gt. 0) then
+       if (.not.(dt_controller%if_variable_dt) .or. &
+       (dt_controller%dt_last_change .gt. this%activ_step - 1)) then
+          call this%project_back(x_u, x_v, x_w, Ax, coef, bclst_u, bclst_v, bclst_w, gs_h, n)
+       end if
+    end if
+
+  end subroutine projection_post_solving_vel
+
+  subroutine bcknd_project_back_vel(this, x_u, x_v, x_w, &
+                                    Ax, coef, bclst_u, bclst_v, bclst_w, gs_h, n)
+    class(projection_vel_t) :: this
+    integer, intent(inout) :: n
+    class(Ax_t), intent(inout) :: Ax
+    class(coef_t), intent(inout) :: coef
+    class(bc_list_t), intent(inout) :: bclst_u, bclst_v, bclst_w
+    type(gs_t), intent(inout) :: gs_h
+    real(kind=rp), intent(inout), dimension(n) :: x_u, x_v, x_w
+    type(c_ptr) :: x_u_d, x_v_d, x_w_d
+
+    call profiler_start_region('Project back', 17)
+
+    if (NEKO_BCKND_DEVICE .eq. 1) then
+       x_u_d = device_get_ptr(x_u)
+       x_v_d = device_get_ptr(x_v)
+       x_w_d = device_get_ptr(x_w)
+       if (this%proj_u%m .gt. 0) then ! Restore desired solution
+          call device_add2(x_u_d, this%proj_u%xbar_d, n)
+       end if
+       if (this%proj_v%m .gt. 0) then ! Restore desired solution
+          call device_add2(x_v_d, this%proj_v%xbar_d, n)
+       end if 
+       if (this%proj_w%m .gt. 0) then ! Restore desired solution
+          call device_add2(x_w_d, this%proj_w%xbar_d, n)
+       end if
+
+       if (this%proj_u%m .eq. this%proj_u%L) then
+          this%proj_u%m = 1
+       else
+          this%proj_u%m = min(this%proj_u%m + 1, this%proj_u%L)
+       end if
+       if (this%proj_v%m .eq. this%proj_v%L) then
+          this%proj_v%m = 1
+       else
+          this%proj_v%m = min(this%proj_v%m + 1, this%proj_v%L)
+       end if
+       if (this%proj_w%m .eq. this%proj_w%L) then
+          this%proj_w%m = 1
+       else
+          this%proj_w%m = min(this%proj_w%m + 1, this%proj_w%L)
+       end if
+
+       call device_copy(this%proj_u%xx_d(this%proj_u%m), &
+                        x_u_d,n) ! Update (X,B)
+       call device_copy(this%proj_v%xx_d(this%proj_u%m), &
+                        x_v_d,n) ! Update (X,B)
+       call device_copy(this%proj_w%xx_d(this%proj_u%m), &
+                        x_w_d,n) ! Update (X,B)
+
+    else
+       if (this%proj_u%m .gt. 0) then
+          call add2(x_u, this%proj_u%xbar, n) ! Restore desired solution
+       end if
+       if (this%proj_v%m .gt. 0) then
+          call add2(x_v, this%proj_v%xbar, n) ! Restore desired solution
+       end if
+       if (this%proj_w%m .gt. 0) then
+          call add2(x_w, this%proj_w%xbar, n) ! Restore desired solution
+       end if
+
+       if (this%proj_u%m .eq. this%proj_u%L) then
+          this%proj_u%m = 1
+       else
+          this%proj_u%m = min(this%proj_u%m + 1, this%proj_u%L)
+       end if
+       if (this%proj_v%m .eq. this%proj_v%L) then
+          this%proj_v%m = 1
+       else
+          this%proj_v%m = min(this%proj_v%m + 1, this%proj_v%L)
+       end if
+       if (this%proj_w%m .eq. this%proj_w%L) then
+          this%proj_w%m = 1
+       else
+          this%proj_w%m = min(this%proj_w%m + 1, this%proj_w%L)
+       end if
+
+       call copy(this%proj_u%xx(1, this%proj_u%m), x_u, n)   ! Update (X,B)
+       call copy(this%proj_v%xx(1, this%proj_v%m), x_v, n)   ! Update (X,B)
+       call copy(this%proj_w%xx(1, this%proj_w%m), x_w, n)   ! Update (X,B)
+    end if
+
+    call Ax%compute_vector(this%proj_u%bb(1,this%proj_u%m), &
+                           this%proj_v%bb(1,this%proj_v%m), &
+                           this%proj_w%bb(1,this%proj_w%m), x_u, x_v, x_w, &
+                           coef, coef%msh, coef%Xh)
+
+    call gs_h%gs_op_vector(this%proj_u%bb(1,this%proj_u%m), n, GS_OP_ADD)
+    call gs_h%gs_op_vector(this%proj_v%bb(1,this%proj_v%m), n, GS_OP_ADD)
+    call gs_h%gs_op_vector(this%proj_w%bb(1,this%proj_w%m), n, GS_OP_ADD)
+
+    call bc_list_apply_scalar(bclst_u, this%proj_u%bb(1,this%proj_u%m), n)
+    call bc_list_apply_scalar(bclst_v, this%proj_v%bb(1,this%proj_v%m), n)
+    call bc_list_apply_scalar(bclst_w, this%proj_w%bb(1,this%proj_w%m), n)
+
+    if (NEKO_BCKND_DEVICE .eq. 1)  then
+       call device_proj_ortho(this%proj_u, this%proj_u%xx_d, &
+                              this%proj_u%bb_d, coef%mult_d, n)
+       call device_proj_ortho(this%proj_v, this%proj_v%xx_d, &
+                              this%proj_v%bb_d, coef%mult_d, n)
+       call device_proj_ortho(this%proj_w, this%proj_w%xx_d, &
+                              this%proj_w%bb_d, coef%mult_d, n)
+    else
+       call cpu_proj_ortho(this%proj_u, this%proj_u%xx, &
+                           this%proj_u%bb, coef%mult, n)
+       call cpu_proj_ortho(this%proj_v, this%proj_v%xx, &
+                           this%proj_v%bb, coef%mult, n)
+       call cpu_proj_ortho(this%proj_w, this%proj_w%xx, &
+                           this%proj_w%bb, coef%mult, n)
+    end if
+    call profiler_end_region('Project back', 17)
+  end subroutine bcknd_project_back_vel
+end module projection_vel

src/fluid/fluid_pnpn.f90

@@ -46,6 +46,7 @@ module fluid_pnpn
   use fluid_aux, only : fluid_step_info
   use time_scheme_controller, only : time_scheme_controller_t
   use projection, only : projection_t
+  use projection_vel, only : projection_vel_t
   use device, only : device_memcpy, HOST_TO_DEVICE
   use advection, only : advection_t, advection_factory
   use profiler, only : profiler_start_region, profiler_end_region
@@ -84,9 +85,7 @@ module fluid_pnpn
      class(ax_t), allocatable :: Ax_prs
 
      type(projection_t) :: proj_prs
-     type(projection_t) :: proj_u
-     type(projection_t) :: proj_v
-     type(projection_t) :: proj_w
+     type(projection_vel_t) :: proj_vel
 
      type(facet_normal_t) :: bc_prs_surface !< Surface term in pressure rhs
      type(facet_normal_t) :: bc_sym_surface !< Surface term in pressure rhs
@@ -345,21 +344,17 @@ subroutine fluid_pnpn_init(this, msh, lx, params, user, time_scheme)
 
     !Intialize projection space thingy
 
-    if (this%variable_material_properties .and. &
-          this%vel_projection_dim .gt. 0) then
-       call neko_error("Velocity projection not available for full stress &
-             &formulation")
-    end if
+   !  if (this%variable_material_properties .and. &
+   !        this%vel_projection_dim .gt. 0) then
+   !     call neko_error("Velocity projection not available for full stress &
+   !           &formulation")
+   !  end if
 
 
     call this%proj_prs%init(this%dm_Xh%size(), this%pr_projection_dim, &
                               this%pr_projection_activ_step)
 
-    call this%proj_u%init(this%dm_Xh%size(), this%vel_projection_dim, &
-                              this%vel_projection_activ_step)
-    call this%proj_v%init(this%dm_Xh%size(), this%vel_projection_dim, &
-                              this%vel_projection_activ_step)
-    call this%proj_w%init(this%dm_Xh%size(), this%vel_projection_dim, &
+    call this%proj_vel%init(this%dm_Xh%size(), this%vel_projection_dim, &
                               this%vel_projection_activ_step)
 
 
@@ -534,9 +529,7 @@ subroutine fluid_pnpn_free(this)
     call bc_list_free(this%bclst_vel_res)
     call bc_list_free(this%bclst_dp)
     call this%proj_prs%free()
-    call this%proj_u%free()
-    call this%proj_v%free()
-    call this%proj_w%free()
+    call this%proj_vel%free()
 
     call this%p_res%free()
     call this%u_res%free()
@@ -786,9 +779,8 @@ subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
 
       call profiler_end_region('Velocity_residual', 19)
 
-      call this%proj_u%pre_solving(u_res%x, tstep, c_Xh, n, dt_controller)
-      call this%proj_v%pre_solving(v_res%x, tstep, c_Xh, n, dt_controller)
-      call this%proj_w%pre_solving(w_res%x, tstep, c_Xh, n, dt_controller)
+      call this%proj_vel%pre_solving(u_res%x, v_res%x, w_res%x, &
+                                     tstep, c_Xh, n, dt_controller)
 
       call this%pc_vel%update()
 
@@ -799,12 +791,8 @@ subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
            this%ksp_vel%max_iter)
       call profiler_end_region("Velocity_solve", 4)
 
-      call this%proj_u%post_solving(du%x, Ax_vel, c_Xh, &
-                                 this%bclst_du, gs_Xh, n, tstep, dt_controller)
-      call this%proj_v%post_solving(dv%x, Ax_vel, c_Xh, &
-                                 this%bclst_dv, gs_Xh, n, tstep, dt_controller)
-      call this%proj_w%post_solving(dw%x, Ax_vel, c_Xh, &
-                                 this%bclst_dw, gs_Xh, n, tstep, dt_controller)
+      call this%proj_vel%post_solving(du%x, dv%x, dw%x, Ax_vel, c_Xh, &
+                                 this%bclst_du, this%bclst_dv, this%bclst_dw, gs_Xh, n, tstep, dt_controller)
 
       if (NEKO_BCKND_DEVICE .eq. 1) then
          call device_opadd2cm(u%x_d, v%x_d, w%x_d, &