Loading [MathJax]/extensions/tex2jax.js
Neko 0.9.99
A portable framework for high-order spectral element flow simulations
All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Pages
projection_vel.f90
Go to the documentation of this file.
1! Copyright (c) 2024, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
36 use num_types, only : rp, c_rp
37 use math, only : add2, copy
38 use coefs, only : coef_t
39 use ax_product, only : ax_t
40 use bc_list, only : bc_list_t
41 use gather_scatter, only : gs_t, gs_op_add
43 use device, only : device_get_ptr
46 use, intrinsic :: iso_c_binding
49
50 implicit none
51 private
52
53 type, public :: projection_vel_t
54 type(projection_t) :: proj_u, proj_v, proj_w
55 integer :: activ_step
56 integer :: l
57 contains
58 procedure, pass(this) :: init => projection_init_vel
59 procedure, pass(this) :: free => projection_free_vel
60 procedure, pass(this) :: pre_solving => projection_pre_solving_vel
61 procedure, pass(this) :: post_solving => projection_post_solving_vel
62 procedure, pass(this) :: project_back => bcknd_project_back_vel
63 end type projection_vel_t
64
65contains
66
67 subroutine projection_init_vel(this, n, L, activ_step)
68 class(projection_vel_t), target, intent(inout) :: this
69 integer, intent(in) :: n
70 integer, optional, intent(in) :: L, activ_step
71 integer :: i
72 integer(c_size_t) :: ptr_size
73 type(c_ptr) :: ptr
74 real(c_rp) :: dummy
75
76 call this%free()
77
78 call this%proj_u%init(n, l, activ_step)
79 call this%proj_v%init(n, l, activ_step)
80 call this%proj_w%init(n, l, activ_step)
81
82 this%L = l
83 this%activ_step = activ_step
84
85
86 end subroutine projection_init_vel
87
88 subroutine projection_free_vel(this)
89 class(projection_vel_t), intent(inout) :: this
90 integer :: i
91
92 call this%proj_u%free()
93 call this%proj_v%free()
94 call this%proj_w%free()
95
96 end subroutine projection_free_vel
97
98 subroutine projection_pre_solving_vel(this, b_u, b_v, b_w, tstep, coef, n, &
99 dt_controller, &
100 stringx, stringy, stringz)
101 class(projection_vel_t), intent(inout) :: this
102 integer, intent(inout) :: n
103 real(kind=rp), intent(inout), dimension(n) :: b_u, b_v, b_w
104 integer, intent(in) :: tstep
105 class(coef_t), intent(inout) :: coef
106 type(time_step_controller_t), intent(in) :: dt_controller
107 character(len=*), optional :: stringx, stringy, stringz
108
109 call this%proj_u%pre_solving(b_u, tstep, coef, n, dt_controller, stringx)
110 call this%proj_v%pre_solving(b_v, tstep, coef, n, dt_controller, stringy)
111 call this%proj_w%pre_solving(b_w, tstep, coef, n, dt_controller, stringz)
112
113 end subroutine projection_pre_solving_vel
114
115 subroutine projection_post_solving_vel(this, x_u, x_v, x_w, Ax, coef, &
116 bclst_u, bclst_v, bclst_w, &
117 gs_h, n, tstep, dt_controller)
118 class(projection_vel_t), intent(inout) :: this
119 integer, intent(inout) :: n
120 class(ax_t), intent(inout) :: Ax
121 class(coef_t), intent(inout) :: coef
122 class(bc_list_t), intent(inout) :: bclst_u, bclst_v, bclst_w
123 type(gs_t), intent(inout) :: gs_h
124 real(kind=rp), intent(inout), dimension(n) :: x_u, x_v, x_w
125 integer, intent(in) :: tstep
126 type(time_step_controller_t), intent(in) :: dt_controller
127
128 ! Here we assume the projection space sizes and activate steps
129 ! for all three velocity equations are the same
130 if (tstep .gt. this%activ_step .and. this%L .gt. 0) then
131 if (.not.(dt_controller%if_variable_dt) .or. &
132 (dt_controller%dt_last_change .gt. this%activ_step - 1)) then
133 call this%project_back(x_u, x_v, x_w, ax, coef, bclst_u, bclst_v, bclst_w, gs_h, n)
134 end if
135 end if
136
137 end subroutine projection_post_solving_vel
138
139 subroutine bcknd_project_back_vel(this, x_u, x_v, x_w, &
140 Ax, coef, bclst_u, bclst_v, bclst_w, gs_h, n)
141 class(projection_vel_t) :: this
142 integer, intent(inout) :: n
143 class(ax_t), intent(inout) :: Ax
144 class(coef_t), intent(inout) :: coef
145 class(bc_list_t), intent(inout) :: bclst_u, bclst_v, bclst_w
146 type(gs_t), intent(inout) :: gs_h
147 real(kind=rp), intent(inout), dimension(n) :: x_u, x_v, x_w
148 type(c_ptr) :: x_u_d, x_v_d, x_w_d
149
150 call profiler_start_region('Project back', 17)
151
152 if (neko_bcknd_device .eq. 1) then
153 x_u_d = device_get_ptr(x_u)
154 x_v_d = device_get_ptr(x_v)
155 x_w_d = device_get_ptr(x_w)
156 if (this%proj_u%m .gt. 0) then ! Restore desired solution
157 call device_add2(x_u_d, this%proj_u%xbar_d, n)
158 end if
159 if (this%proj_v%m .gt. 0) then ! Restore desired solution
160 call device_add2(x_v_d, this%proj_v%xbar_d, n)
161 end if
162 if (this%proj_w%m .gt. 0) then ! Restore desired solution
163 call device_add2(x_w_d, this%proj_w%xbar_d, n)
164 end if
165
166 if (this%proj_u%m .eq. this%proj_u%L) then
167 this%proj_u%m = 1
168 else
169 this%proj_u%m = min(this%proj_u%m + 1, this%proj_u%L)
170 end if
171 if (this%proj_v%m .eq. this%proj_v%L) then
172 this%proj_v%m = 1
173 else
174 this%proj_v%m = min(this%proj_v%m + 1, this%proj_v%L)
175 end if
176 if (this%proj_w%m .eq. this%proj_w%L) then
177 this%proj_w%m = 1
178 else
179 this%proj_w%m = min(this%proj_w%m + 1, this%proj_w%L)
180 end if
181
182 call device_copy(this%proj_u%xx_d(this%proj_u%m), &
183 x_u_d,n) ! Update (X,B)
184 call device_copy(this%proj_v%xx_d(this%proj_u%m), &
185 x_v_d,n) ! Update (X,B)
186 call device_copy(this%proj_w%xx_d(this%proj_u%m), &
187 x_w_d,n) ! Update (X,B)
188
189 else
190 if (this%proj_u%m .gt. 0) then
191 call add2(x_u, this%proj_u%xbar, n) ! Restore desired solution
192 end if
193 if (this%proj_v%m .gt. 0) then
194 call add2(x_v, this%proj_v%xbar, n) ! Restore desired solution
195 end if
196 if (this%proj_w%m .gt. 0) then
197 call add2(x_w, this%proj_w%xbar, n) ! Restore desired solution
198 end if
199
200 if (this%proj_u%m .eq. this%proj_u%L) then
201 this%proj_u%m = 1
202 else
203 this%proj_u%m = min(this%proj_u%m + 1, this%proj_u%L)
204 end if
205 if (this%proj_v%m .eq. this%proj_v%L) then
206 this%proj_v%m = 1
207 else
208 this%proj_v%m = min(this%proj_v%m + 1, this%proj_v%L)
209 end if
210 if (this%proj_w%m .eq. this%proj_w%L) then
211 this%proj_w%m = 1
212 else
213 this%proj_w%m = min(this%proj_w%m + 1, this%proj_w%L)
214 end if
215
216 call copy(this%proj_u%xx(1, this%proj_u%m), x_u, n) ! Update (X,B)
217 call copy(this%proj_v%xx(1, this%proj_v%m), x_v, n) ! Update (X,B)
218 call copy(this%proj_w%xx(1, this%proj_w%m), x_w, n) ! Update (X,B)
219 end if
220
221 call ax%compute_vector(this%proj_u%bb(1,this%proj_u%m), &
222 this%proj_v%bb(1,this%proj_v%m), &
223 this%proj_w%bb(1,this%proj_w%m), x_u, x_v, x_w, &
224 coef, coef%msh, coef%Xh)
225
226 call gs_h%gs_op_vector(this%proj_u%bb(1,this%proj_u%m), n, gs_op_add)
227 call gs_h%gs_op_vector(this%proj_v%bb(1,this%proj_v%m), n, gs_op_add)
228 call gs_h%gs_op_vector(this%proj_w%bb(1,this%proj_w%m), n, gs_op_add)
229
230 call bclst_u%apply_scalar(this%proj_u%bb(1,this%proj_u%m), n)
231 call bclst_v%apply_scalar(this%proj_v%bb(1,this%proj_v%m), n)
232 call bclst_w%apply_scalar(this%proj_w%bb(1,this%proj_w%m), n)
233
234 call proj_ortho(this%proj_u, coef, n)
235 call proj_ortho(this%proj_v, coef, n)
236 call proj_ortho(this%proj_w, coef, n)
237 call profiler_end_region('Project back', 17)
238 end subroutine bcknd_project_back_vel
239end module projection_vel
Return the device pointer for an associated Fortran array.
Definition device.F90:95
Defines a Matrix-vector product.
Definition ax.f90:34
Defines a list of bc_t.
Definition bc_list.f90:34
Coefficients.
Definition coef.f90:34
subroutine, public device_add2(a_d, b_d, n)
Vector addition .
subroutine, public device_copy(a_d, b_d, n)
Copy a vector .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
Gather-scatter.
Definition math.f90:60
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:586
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:238
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public c_rp
Definition num_types.f90:13
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Profiling interface.
Definition profiler.F90:34
subroutine, public profiler_start_region(name, region_id)
Started a named (name) profiler region.
Definition profiler.F90:78
subroutine, public profiler_end_region(name, region_id)
End the most recently started profiler region.
Definition profiler.F90:109
Project x onto X , the space of old solutions and back again Couple projections for velocity.
subroutine projection_init_vel(this, n, l, activ_step)
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)
subroutine bcknd_project_back_vel(this, x_u, x_v, x_w, ax, coef, bclst_u, bclst_v, bclst_w, gs_h, n)
subroutine projection_free_vel(this)
subroutine projection_pre_solving_vel(this, b_u, b_v, b_w, tstep, coef, n, dt_controller, stringx, stringy, stringz)
Project x onto X, the space of old solutions and back again.
subroutine, public proj_ortho(this, coef, n)
Implements type time_step_controller.
Base type for a matrix-vector product providing .
Definition ax.f90:43
A list of allocatable `bc_t`. Follows the standard interface of lists.
Definition bc_list.f90:47
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55