Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
PDE_filter.f90
Go to the documentation of this file.
1! Copyright (c) 2023, 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!
33!
35
37 use num_types, only: rp
38 use json_module, only: json_file
41 use field, only: field_t
42 use coefs, only: coef_t
43 use ax_product, only: ax_t, ax_helm_factory
44 use krylov, only: ksp_t, ksp_monitor_t, krylov_solver_factory
45 use precon, only: pc_t, precon_factory, precon_destroy
46 use bc_list, only : bc_list_t
47 use neumann, only: neumann_t
49 use gather_scatter, only: gs_t, gs_op_add
53 use filter, only: filter_t
56 use coefs, only: coef_t
57 use logger, only: neko_log, log_size
59 use dofmap, only: dofmap_t
60 use jacobi, only: jacobi_t
62 use sx_jacobi, only: sx_jacobi_t
63 use utils, only: neko_error
65 implicit none
66 private
67
72 type, public, extends(filter_t) :: pde_filter_t
73
75 class(ax_t), allocatable :: ax
77 type(ksp_monitor_t) :: ksp_results(1)
79 class(ksp_t), allocatable :: ksp_filt
81 class(pc_t), allocatable :: pc_filt
83 type(bc_list_t) :: bclst_filt
84
85 ! Inputs from the user
87 real(kind=rp) :: r
89 real(kind=rp) :: abstol_filt
91 integer :: ksp_max_iter
93 character(len=:), allocatable :: ksp_solver
94 ! > preconditioner type
95 character(len=:), allocatable :: precon_type_filt
96 integer :: ksp_n, n, i
97
98
99
100 contains
102 procedure, pass(this) :: init => pde_filter_init_from_json
104 procedure, pass(this) :: init_from_components => &
107 procedure, pass(this) :: free => pde_filter_free
109 procedure, pass(this) :: apply => pde_filter_apply
110 end type pde_filter_t
111
112contains
113
115 subroutine pde_filter_init_from_json(this, json, coef)
116 class(pde_filter_t), intent(inout) :: this
117 type(json_file), intent(inout) :: json
118 type(coef_t), intent(in) :: coef
119
120 ! user parameters
121 call json_get(json, "filter.radius", this%r)
122
123 call json_get_or_default(json, "filter.tolerance", this%abstol_filt, &
124 1.0e-10_rp)
125
126 call json_get_or_default(json, "filter.max_iter", this%ksp_max_iter, 200)
127
128 call json_get_or_default(json, "filter.solver", this%ksp_solver, 'cg')
129
130 call json_get_or_default(json, "filter.preconditioner", &
131 this%precon_type_filt, 'jacobi')
132
133 call this%init_base(json, coef)
134 call pde_filter_init_from_components(this, coef)
135
136 end subroutine pde_filter_init_from_json
137
139 subroutine pde_filter_init_from_components(this, coef)
140 class(pde_filter_t), intent(inout) :: this
141 type(coef_t), intent(in) :: coef
142 integer :: n
143
144 n = this%coef%dof%size()
145
146 ! init the bc list (all Neuman BCs, will remain empty)
147 call this%bclst_filt%init()
148
149 ! Setup backend dependent Ax routines
150 call ax_helm_factory(this%Ax, full_formulation = .false.)
151
152 ! set up krylov solver
153 call krylov_solver_factory(this%ksp_filt, n, this%ksp_solver, &
154 this%ksp_max_iter, this%abstol_filt)
155
156 ! set up preconditioner
157 call filter_precon_factory(this%pc_filt, this%ksp_filt, &
158 this%coef, this%coef%dof, &
159 this%coef%gs_h, &
160 this%bclst_filt, this%precon_type_filt)
161
163
165 subroutine pde_filter_free(this)
166 class(pde_filter_t), intent(inout) :: this
167
168 if (allocated(this%Ax)) then
169 deallocate(this%Ax)
170 end if
171
172 if (allocated(this%ksp_filt)) then
173 call this%ksp_filt%free()
174 deallocate(this%ksp_filt)
175 end if
176
177 if (allocated(this%pc_filt)) then
178 call precon_destroy(this%pc_filt)
179 deallocate(this%pc_filt)
180 end if
181
182 call this%bclst_filt%free()
183
184 call this%free_base()
185
186 end subroutine pde_filter_free
187
191 subroutine pde_filter_apply(this, F_out, F_in)
192 class(pde_filter_t), intent(inout) :: this
193 type(field_t), intent(in) :: F_in
194 type(field_t), intent(inout) :: F_out
195 integer :: n, i
196 ! type(field_t), pointer :: RHS
197 type(field_t) :: RHS, d_F_out
198 character(len=LOG_SIZE) :: log_buf
199 ! integer :: temp_indices(1)
200
201 n = this%coef%dof%size()
202 ! TODO
203 ! This is a bit awkward, because the init for the source terms occurs
204 ! before the init of the scratch registry.
205 ! So we can't use the scratch registry here.
206 ! call neko_scratch_registry%request_field(RHS, temp_indices(1))
207 call rhs%init(this%coef%dof)
208 call d_f_out%init(this%coef%dof)
209
210 ! in a similar fasion to pressure/velocity, we will solve for d_F_out.
211
212 ! to improve convergence, we use F_in as an initial guess for F_out.
213 ! so F_out = F_in + d_F_in.
214
215 ! Defining the operator A = -r^2 \nabla^2 + I
216 ! the system changes from:
217 ! A (F_out) = F_in
218 ! to
219 ! A (d_F_out) = F_in - A(F_in)
220
221 ! set up Helmholtz operators and RHS
222 if (neko_bcknd_device .eq. 1) then
223 ! TODO
224 ! I think this is correct but I've never tested it
225 call device_cfill(this%coef%h1_d, this%r**2, n)
226 call device_cfill(this%coef%h2_d, 1.0_rp, n)
227 else
228 do i = 1, n
229 ! h1 is already negative in its definition
230 this%coef%h1(i,1,1,1) = this%r**2
231 ! ax_helm includes the mass matrix in h2
232 this%coef%h2(i,1,1,1) = 1.0_rp
233 end do
234 end if
235 this%coef%ifh2 = .true.
236
237 ! compute the A(F_in) component of the RHS
238 ! (note, to be safe with the inout intent we first copy F_in to the
239 ! temporary d_F_out)
240 call field_copy(d_f_out, f_in)
241 call this%Ax%compute(rhs%x, d_f_out%x, this%coef, this%coef%msh, &
242 this%coef%Xh)
243
244 if (neko_bcknd_device .eq. 1) then
245 call device_subcol3(rhs%x_d, f_in%x_d, this%coef%B_d, n)
246 call device_cmult(rhs%x_d, -1.0_rp, n)
247 else
248 do i = 1, n
249 ! mass matrix should be included here
250 rhs%x(i,1,1,1) = f_in%x(i,1,1,1) * this%coef%B(i,1,1,1) &
251 - rhs%x(i,1,1,1)
252 end do
253 end if
254
255 ! gather scatter
256 call this%coef%gs_h%op(rhs, gs_op_add)
257
258 ! set BCs
259 call this%bclst_filt%apply_scalar(rhs%x, n)
260
261 ! Solve Helmholtz equation
262 call profiler_start_region('filter solve')
263 this%ksp_results(1) = &
264 this%ksp_filt%solve(this%Ax, d_f_out, rhs%x, n, this%coef, &
265 this%bclst_filt, this%coef%gs_h)
266
268
269 ! add result
270 call field_add3(f_out, f_in, d_f_out)
271 ! update preconditioner (needed?)
272 call this%pc_filt%update()
273
274 ! write it all out
275 call neko_log%message('Filter')
276
277 write(log_buf, '(A,A,A)') 'Iterations: ',&
278 'Start residual: ', 'Final residual:'
279 call neko_log%message(log_buf)
280 write(log_buf, '(I11,3x, E15.7,5x, E15.7)') this%ksp_results%iter, &
281 this%ksp_results%res_start, this%ksp_results%res_final
282 call neko_log%message(log_buf)
283
284 !call neko_scratch_registry%relinquish_field(temp_indices)
285 call rhs%free()
286 call d_f_out%free()
287
288 end subroutine pde_filter_apply
289
291 subroutine filter_precon_factory(pc, ksp, coef, dof, gs, bclst, &
292 pctype)
293 class(pc_t), allocatable, target, intent(inout) :: pc
294 class(ksp_t), target, intent(inout) :: ksp
295 type(coef_t), target, intent(in) :: coef
296 type(dofmap_t), target, intent(in) :: dof
297 type(gs_t), target, intent(inout) :: gs
298 type(bc_list_t), target, intent(inout) :: bclst
299 character(len=*) :: pctype
300
301 call precon_factory(pc, pctype)
302
303 select type (pcp => pc)
304 type is (jacobi_t)
305 call pcp%init(coef, dof, gs)
306 type is (sx_jacobi_t)
307 call pcp%init(coef, dof, gs)
308 type is (device_jacobi_t)
309 call pcp%init(coef, dof, gs)
310 end select
311
312 call ksp%set_pc(pc)
313
314 end subroutine filter_precon_factory
315
316end module pde_filter
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Retrieves a parameter by name or throws an error.
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
Jacobi preconditioner accelerator backend.
subroutine, public device_cmult(a_d, c, n, strm)
Multiplication by constant c .
subroutine, public device_subcol3(a_d, b_d, c_d, n, strm)
Returns .
subroutine, public device_cfill(a_d, c, n, strm)
Set all elements to a constant c .
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
subroutine, public field_copy(a, b, n)
Copy a vector .
subroutine, public field_add3(a, b, c, n)
Vector addition .
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Defines a field.
Definition field.f90:34
Filter to be applied to a scalar field.
Definition filter.f90:38
Gather-scatter.
Jacobi preconditioner.
Definition pc_jacobi.f90:34
Utilities for retrieving parameters from the case files.
Implements the base abstract type for Krylov solvers plus helper types.
Definition krylov.f90:34
integer, parameter, public ksp_max_iter
Maximum number of iters.
Definition krylov.f90:51
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
Defines a mesh.
Definition mesh.f90:34
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
Definition mesh.f90:62
integer, parameter, public neko_msh_max_zlbl_len
Max length of a zone label.
Definition mesh.f90:64
Build configurations.
integer, parameter neko_bcknd_device
Defines a Neumann boundary condition.
Definition neumann.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
A PDE based filter.
subroutine pde_filter_free(this)
Destructor.
subroutine pde_filter_init_from_json(this, json, coef)
Constructor from json.
subroutine filter_precon_factory(pc, ksp, coef, dof, gs, bclst, pctype)
Initialize a Krylov preconditioner.
subroutine pde_filter_init_from_components(this, coef)
Actual constructor.
subroutine pde_filter_apply(this, f_out, f_in)
Apply the filter.
Defines Pressure and velocity residuals in the Pn-Pn formulation.
Definition pnpn_res.f90:34
Krylov preconditioner.
Definition precon.f90:34
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
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
type(scratch_registry_t), target, public neko_scratch_registry
Global scratch registry.
Jacobi preconditioner SX-Aurora backend.
Utilities.
Definition utils.f90:35
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:48
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Defines a jacobi preconditioner.
Base abstract class for filter.
Definition filter.f90:48
Defines a jacobi preconditioner.
Definition pc_jacobi.f90:45
Type for storing initial and final residuals in a Krylov solver.
Definition krylov.f90:56
Base abstract type for a canonical Krylov method, solving .
Definition krylov.f90:73
A Neumann boundary condition for scalar fields. Sets the flux of the field to the chosen value.
Definition neumann.f90:56
A PDE based filter mapping $\rho \mapsto \tilde{\rho}$, see Lazarov & O. Sigmund 2010,...
Abstract type to compute pressure residual.
Definition pnpn_res.f90:48
Defines a canonical Krylov preconditioner.
Definition precon.f90:40
Defines a jacobi preconditioner for SX-Aurora.