Neko  0.9.0
A portable framework for high-order spectral element flow simulations
fluid_volflow.f90
Go to the documentation of this file.
1 ! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
2 !
3 ! The UChicago Argonne, LLC as Operator of Argonne National
4 ! Laboratory holds copyright in the Software. The copyright holder
5 ! reserves all rights except those expressly granted to licensees,
6 ! and U.S. Government license rights.
7 !
8 ! Redistribution and use in source and binary forms, with or without
9 ! modification, are permitted provided that the following conditions
10 ! are met:
11 !
12 ! 1. Redistributions of source code must retain the above copyright
13 ! notice, this list of conditions and the disclaimer below.
14 !
15 ! 2. Redistributions in binary form must reproduce the above copyright
16 ! notice, this list of conditions and the disclaimer (as noted below)
17 ! in the documentation and/or other materials provided with the
18 ! distribution.
19 !
20 ! 3. Neither the name of ANL nor the names of its contributors
21 ! may be used to endorse or promote products derived from this software
22 ! without specific prior written permission.
23 !
24 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
27 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
28 ! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
29 ! ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30 ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
31 ! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32 ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33 ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34 ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35 ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 !
37 ! Additional BSD Notice
38 ! ---------------------
39 ! 1. This notice is required to be provided under our contract with
40 ! the U.S. Department of Energy (DOE). This work was produced at
41 ! Argonne National Laboratory under Contract
42 ! No. DE-AC02-06CH11357 with the DOE.
43 !
44 ! 2. Neither the United States Government nor UCHICAGO ARGONNE,
45 ! LLC nor any of their employees, makes any warranty,
46 ! express or implied, or assumes any liability or responsibility for the
47 ! accuracy, completeness, or usefulness of any information, apparatus,
48 ! product, or process disclosed, or represents that its use would not
49 ! infringe privately-owned rights.
50 !
51 ! 3. Also, reference herein to any specific commercial products, process,
52 ! or services by trade name, trademark, manufacturer or otherwise does
53 ! not necessarily constitute or imply its endorsement, recommendation,
54 ! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
55 ! The views and opinions of authors expressed
56 ! herein do not necessarily state or reflect those of the United States
57 ! Government or UCHICAGO ARGONNE, LLC, and shall
58 ! not be used for advertising or product endorsement purposes.
59 !
61  use operators, only : opgrad, cdtp
62  use num_types, only : rp
63  use mathops, only : opchsign
64  use krylov, only : ksp_t, ksp_monitor_t
65  use precon, only : pc_t
66  use dofmap, only : dofmap_t
67  use field, only : field_t
68  use coefs, only : coef_t
70  use math, only : copy, glsc2, glmin, glmax, add2, add2s2
71  use comm
72  use neko_config, only : neko_bcknd_device
76  use gather_scatter, only : gs_t, gs_op_add
77  use json_module, only : json_file
78  use json_utils, only: json_get
81  use ax_product, only : ax_t
82  implicit none
83  private
84 
86  type, public :: fluid_volflow_t
87  integer :: flow_dir
88  logical :: avflow
89  real(kind=rp) :: flow_rate
90  real(kind=rp) :: dtlag = 0d0
91  real(kind=rp) :: bdlag = 0d0
92  type(field_t) :: u_vol, v_vol, w_vol, p_vol
93  real(kind=rp) :: domain_length, base_flow
95  type(scratch_registry_t) :: scratch
96  contains
97  procedure, pass(this) :: init => fluid_vol_flow_init
98  procedure, pass(this) :: free => fluid_vol_flow_free
99  procedure, pass(this) :: adjust => fluid_vol_flow
100  procedure, private, pass(this) :: compute => fluid_vol_flow_compute
101  end type fluid_volflow_t
102 
103 contains
104 
105  subroutine fluid_vol_flow_init(this, dm_Xh, params)
106  class(fluid_volflow_t), intent(inout) :: this
107  type(dofmap_t), target, intent(inout) :: dm_Xh
108  type(json_file), intent(inout) :: params
109  logical average
110  integer :: direction
111  real(kind=rp) :: rate
112 
113  call this%free()
114 
115  !Initialize vol_flow (if there is a forced volume flow)
116  call json_get(params, 'case.fluid.flow_rate_force.direction', direction)
117  call json_get(params, 'case.fluid.flow_rate_force.value', rate)
118  call json_get(params, 'case.fluid.flow_rate_force.use_averaged_flow',&
119  average)
120 
121  this%flow_dir = direction
122  this%avflow = average
123  this%flow_rate = rate
124 
125  if (this%flow_dir .ne. 0) then
126  call this%u_vol%init(dm_xh, 'u_vol')
127  call this%v_vol%init(dm_xh, 'v_vol')
128  call this%w_vol%init(dm_xh, 'w_vol')
129  call this%p_vol%init(dm_xh, 'p_vol')
130  end if
131 
132  this%scratch = scratch_registry_t(dm_xh, 3, 1)
133 
134  end subroutine fluid_vol_flow_init
135 
136  subroutine fluid_vol_flow_free(this)
137  class(fluid_volflow_t), intent(inout) :: this
138 
139  call this%u_vol%free()
140  call this%v_vol%free()
141  call this%w_vol%free()
142  call this%p_vol%free()
143 
144  call this%scratch%free()
145 
146  end subroutine fluid_vol_flow_free
147 
151  subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, &
152  ext_bdf, gs_Xh, c_Xh, rho, mu, bd, dt, &
153  bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
154  Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
155  class(fluid_volflow_t), intent(inout) :: this
156  type(field_t), intent(inout) :: u_res, v_res, w_res, p_res
157  type(coef_t), intent(inout) :: c_Xh
158  type(gs_t), intent(inout) :: gs_Xh
159  type(time_scheme_controller_t), intent(inout) :: ext_bdf
160  type(bc_list_t), intent(inout) :: bclst_dp, bclst_du, bclst_dv, bclst_dw
161  type(bc_list_t), intent(inout) :: bclst_vel_res
162  class(ax_t), intent(inout) :: Ax_vel
163  class(ax_t), intent(inout) :: Ax_prs
164  class(ksp_t), intent(inout) :: ksp_prs, ksp_vel
165  class(pc_t), intent(inout) :: pc_prs, pc_vel
166  real(kind=rp), intent(inout) :: bd
167  real(kind=rp), intent(in) :: rho, mu, dt
168  integer, intent(in) :: vel_max_iter, prs_max_iter
169  integer :: n, i
170  real(kind=rp) :: xlmin, xlmax
171  real(kind=rp) :: ylmin, ylmax
172  real(kind=rp) :: zlmin, zlmax
173  type(ksp_monitor_t) :: ksp_results(4)
174  type(field_t), pointer :: ta1, ta2, ta3
175  integer :: temp_indices(3)
176 
177  call this%scratch%request_field(ta1, temp_indices(1))
178  call this%scratch%request_field(ta2, temp_indices(2))
179  call this%scratch%request_field(ta3, temp_indices(3))
180 
181 
182  associate(msh => c_xh%msh, p_vol => this%p_vol, &
183  u_vol => this%u_vol, v_vol => this%v_vol, w_vol => this%w_vol)
184 
185  n = c_xh%dof%size()
186  xlmin = glmin(c_xh%dof%x, n)
187  xlmax = glmax(c_xh%dof%x, n)
188  ylmin = glmin(c_xh%dof%y, n) ! for Y!
189  ylmax = glmax(c_xh%dof%y, n)
190  zlmin = glmin(c_xh%dof%z, n) ! for Z!
191  zlmax = glmax(c_xh%dof%z, n)
192  if (this%flow_dir .eq. 1) then
193  this%domain_length = xlmax - xlmin
194  end if
195  if (this%flow_dir .eq. 2) then
196  this%domain_length = ylmax - ylmin
197  end if
198  if (this%flow_dir .eq. 3) then
199  this%domain_length = zlmax - zlmin
200  end if
201 
202  if (neko_bcknd_device .eq. 1) then
203  call device_cfill(c_xh%h1_d, 1.0_rp/rho, n)
204  call device_rzero(c_xh%h2_d, n)
205  else
206  do i = 1, n
207  c_xh%h1(i,1,1,1) = 1.0_rp / rho
208  c_xh%h2(i,1,1,1) = 0.0_rp
209  end do
210  end if
211  c_xh%ifh2 = .false.
212 
213  ! Compute pressure
214 
215  if (this%flow_dir .eq. 1) then
216  call cdtp(p_res%x, c_xh%h1, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
217  end if
218 
219  if (this%flow_dir .eq. 2) then
220  call cdtp(p_res%x, c_xh%h1, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
221  end if
222 
223  if (this%flow_dir .eq. 3) then
224  call cdtp(p_res%x, c_xh%h1, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
225  end if
226 
227  call gs_xh%op(p_res, gs_op_add)
228  call bc_list_apply_scalar(bclst_dp, p_res%x, n)
229  call pc_prs%update()
230  ksp_results(1) = ksp_prs%solve(ax_prs, p_vol, p_res%x, n, &
231  c_xh, bclst_dp, gs_xh, prs_max_iter)
232 
233  ! Compute velocity
234 
235  call opgrad(u_res%x, v_res%x, w_res%x, p_vol%x, c_xh)
236 
237  if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
238  (neko_bcknd_opencl .eq. 1)) then
239  call device_opchsign(u_res%x_d, v_res%x_d, w_res%x_d, msh%gdim, n)
240  call device_copy(ta1%x_d, c_xh%B_d, n)
241  call device_copy(ta2%x_d, c_xh%B_d, n)
242  call device_copy(ta3%x_d, c_xh%B_d, n)
243  else
244  call opchsign(u_res%x, v_res%x, w_res%x, msh%gdim, n)
245  call copy(ta1%x, c_xh%B, n)
246  call copy(ta2%x, c_xh%B, n)
247  call copy(ta3%x, c_xh%B, n)
248  end if
249  call bc_list_apply_vector(bclst_vel_res,&
250  ta1%x, ta2%x, ta3%x, n)
251 
252  ! add forcing
253 
254  if (neko_bcknd_device .eq. 1) then
255  if (this%flow_dir .eq. 1) then
256  call device_add2(u_res%x_d, ta1%x_d, n)
257  else if (this%flow_dir .eq. 2) then
258  call device_add2(v_res%x_d, ta2%x_d, n)
259  else if (this%flow_dir .eq. 3) then
260  call device_add2(w_res%x_d, ta3%x_d, n)
261  end if
262  else
263  if (this%flow_dir .eq. 1) then
264  call add2(u_res%x, ta1%x, n)
265  else if (this%flow_dir .eq. 2) then
266  call add2(v_res%x, ta2%x, n)
267  else if (this%flow_dir .eq. 3) then
268  call add2(w_res%x, ta3%x, n)
269  end if
270  end if
271 
272  if (neko_bcknd_device .eq. 1) then
273  call device_cfill(c_xh%h1_d, mu, n)
274  call device_cfill(c_xh%h2_d, rho * (bd / dt), n)
275  else
276  do i = 1, n
277  c_xh%h1(i,1,1,1) = mu
278  c_xh%h2(i,1,1,1) = rho * (bd / dt)
279  end do
280  end if
281  c_xh%ifh2 = .true.
282 
283  call gs_xh%op(u_res, gs_op_add)
284  call gs_xh%op(v_res, gs_op_add)
285  call gs_xh%op(w_res, gs_op_add)
286 
287  call bc_list_apply_vector(bclst_vel_res,&
288  u_res%x, v_res%x, w_res%x, n)
289  call pc_vel%update()
290 
291  ksp_results(2:4) = ksp_vel%solve_coupled(ax_vel, &
292  u_vol, v_vol, w_vol, &
293  u_res%x, v_res%x, w_res%x, &
294  n, c_xh, &
295  bclst_du, bclst_dv, bclst_dw, &
296  gs_xh, vel_max_iter)
297 
298  if (neko_bcknd_device .eq. 1) then
299  if (this%flow_dir .eq. 1) then
300  this%base_flow = &
301  device_glsc2(u_vol%x_d, c_xh%B_d, n) / this%domain_length
302  end if
303 
304  if (this%flow_dir .eq. 2) then
305  this%base_flow = &
306  device_glsc2(v_vol%x_d, c_xh%B_d, n) / this%domain_length
307  end if
308 
309  if (this%flow_dir .eq. 3) then
310  this%base_flow = &
311  device_glsc2(w_vol%x_d, c_xh%B_d, n) / this%domain_length
312  end if
313  else
314  if (this%flow_dir .eq. 1) then
315  this%base_flow = glsc2(u_vol%x, c_xh%B, n) / this%domain_length
316  end if
317 
318  if (this%flow_dir .eq. 2) then
319  this%base_flow = glsc2(v_vol%x, c_xh%B, n) / this%domain_length
320  end if
321 
322  if (this%flow_dir .eq. 3) then
323  this%base_flow = glsc2(w_vol%x, c_xh%B, n) / this%domain_length
324  end if
325  end if
326  end associate
327 
328  call this%scratch%relinquish_field(temp_indices)
329  end subroutine fluid_vol_flow_compute
330 
340  subroutine fluid_vol_flow(this, u, v, w, p, u_res, v_res, w_res, p_res, &
341  c_Xh, gs_Xh, ext_bdf, rho, mu, dt, &
342  bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
343  Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
344 
345  class(fluid_volflow_t), intent(inout) :: this
346  type(field_t), intent(inout) :: u, v, w, p
347  type(field_t), intent(inout) :: u_res, v_res, w_res, p_res
348  type(coef_t), intent(inout) :: c_Xh
349  type(gs_t), intent(inout) :: gs_Xh
350  type(time_scheme_controller_t), intent(inout) :: ext_bdf
351  real(kind=rp), intent(in) :: rho, mu, dt
352  type(bc_list_t), intent(inout) :: bclst_dp, bclst_du, bclst_dv, bclst_dw
353  type(bc_list_t), intent(inout) :: bclst_vel_res
354  class(ax_t), intent(inout) :: Ax_vel
355  class(ax_t), intent(inout) :: Ax_prs
356  class(ksp_t), intent(inout) :: ksp_prs, ksp_vel
357  class(pc_t), intent(inout) :: pc_prs, pc_vel
358  integer, intent(in) :: prs_max_iter, vel_max_iter
359  real(kind=rp) :: ifcomp, flow_rate, xsec
360  real(kind=rp) :: current_flow, delta_flow, scale
361  integer :: n, ierr, i
362  type(field_t), pointer :: ta1, ta2, ta3
363 
364  associate(u_vol => this%u_vol, v_vol => this%v_vol, &
365  w_vol => this%w_vol, p_vol => this%p_vol)
366 
367  n = c_xh%dof%size()
368 
369  ! If either dt or the backwards difference coefficient change,
370  ! then recompute base flow solution corresponding to unit forcing:
371 
372  ifcomp = 0.0_rp
373 
374  if (dt .ne. this%dtlag .or. &
375  ext_bdf%diffusion_coeffs(1) .ne. this%bdlag) then
376  ifcomp = 1.0_rp
377  end if
378 
379  this%dtlag = dt
380  this%bdlag = ext_bdf%diffusion_coeffs(1)
381 
382  call mpi_allreduce(mpi_in_place, ifcomp, 1, &
383  mpi_real_precision, mpi_sum, neko_comm, ierr)
384 
385  if (ifcomp .gt. 0d0) then
386  call this%compute(u_res, v_res, w_res, p_res, &
387  ext_bdf, gs_xh, c_xh, rho, mu, ext_bdf%diffusion_coeffs(1), dt, &
388  bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, &
389  ax_vel, ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, &
390  vel_max_iter)
391  end if
392 
393  if (neko_bcknd_device .eq. 1) then
394  if (this%flow_dir .eq. 1) then
395  current_flow = &
396  device_glsc2(u%x_d, c_xh%B_d, n) / this%domain_length ! for X
397  else if (this%flow_dir .eq. 2) then
398  current_flow = &
399  device_glsc2(v%x_d, c_xh%B_d, n) / this%domain_length ! for Y
400  else if (this%flow_dir .eq. 3) then
401  current_flow = &
402  device_glsc2(w%x_d, c_xh%B_d, n) / this%domain_length ! for Z
403  end if
404  else
405  if (this%flow_dir .eq. 1) then
406  current_flow = glsc2(u%x, c_xh%B, n) / this%domain_length ! for X
407  else if (this%flow_dir .eq. 2) then
408  current_flow = glsc2(v%x, c_xh%B, n) / this%domain_length ! for Y
409  else if (this%flow_dir .eq. 3) then
410  current_flow = glsc2(w%x, c_xh%B, n) / this%domain_length ! for Z
411  end if
412  end if
413 
414  if (this%avflow) then
415  xsec = c_xh%volume / this%domain_length
416  flow_rate = this%flow_rate*xsec
417  else
418  flow_rate = this%flow_rate
419  end if
420 
421  delta_flow = flow_rate - current_flow
422  scale = delta_flow / this%base_flow
423 
424  if (neko_bcknd_device .eq. 1) then
425  call device_add2s2(u%x_d, u_vol%x_d, scale, n)
426  call device_add2s2(v%x_d, v_vol%x_d, scale, n)
427  call device_add2s2(w%x_d, w_vol%x_d, scale, n)
428  call device_add2s2(p%x_d, p_vol%x_d, scale, n)
429  else
430  do concurrent(i = 1: n)
431  u%x(i,1,1,1) = u%x(i,1,1,1) + scale * u_vol%x(i,1,1,1)
432  v%x(i,1,1,1) = v%x(i,1,1,1) + scale * v_vol%x(i,1,1,1)
433  w%x(i,1,1,1) = w%x(i,1,1,1) + scale * w_vol%x(i,1,1,1)
434  p%x(i,1,1,1) = p%x(i,1,1,1) + scale * p_vol%x(i,1,1,1)
435  end do
436  end if
437  end associate
438 
439  end subroutine fluid_vol_flow
440 
441 
442 end module fluid_volflow
Retrieves a parameter by name or throws an error.
Definition: json_utils.f90:45
Defines a Matrix-vector product.
Definition: ax.f90:34
Defines a boundary condition.
Definition: bc.f90:34
subroutine, public bc_list_apply_vector(bclst, x, y, z, n, t, tstep)
Apply a list of boundary conditions to a vector field.
Definition: bc.f90:587
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
Definition: bc.f90:530
Coefficients.
Definition: coef.f90:34
Definition: comm.F90:1
subroutine, public device_add2(a_d, b_d, n)
Vector addition .
subroutine, public device_rzero(a_d, n)
Zero a real vector.
subroutine, public device_add2s2(a_d, b_d, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
real(kind=rp) function, public device_glsc2(a_d, b_d, n)
Weighted inner product .
subroutine, public device_copy(a_d, b_d, n)
Copy a vector .
Definition: device_math.F90:76
subroutine, public device_cfill(a_d, c, n)
Set all elements to a constant c .
subroutine, public device_opchsign(a1_d, a2_d, a3_d, gdim, n)
Defines a mapping of the degrees of freedom.
Definition: dofmap.f90:35
Defines a field.
Definition: field.f90:34
subroutine fluid_vol_flow_init(this, dm_Xh, params)
subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, ext_bdf, gs_Xh, c_Xh, rho, mu, bd, dt, bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
Compute flow adjustment.
subroutine fluid_vol_flow_free(this)
subroutine fluid_vol_flow(this, u, v, w, p, u_res, v_res, w_res, p_res, c_Xh, gs_Xh, ext_bdf, rho, mu, dt, bclst_dp, bclst_du, bclst_dv, bclst_dw, bclst_vel_res, Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, vel_max_iter)
Adjust flow volume.
Gather-scatter.
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
Implements the base abstract type for Krylov solvers plus helper types.
Definition: krylov.f90:34
Definition: math.f90:60
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
Definition: math.f90:876
subroutine, public add2(a, b, n)
Vector addition .
Definition: math.f90:587
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition: math.f90:377
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
Definition: math.f90:673
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition: math.f90:407
Collection of vector field operations operating on and . Note that in general the indices and ....
Definition: mathops.f90:65
subroutine, public opchsign(a1, a2, a3, gdim, n)
for and .
Definition: mathops.f90:76
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Operators.
Definition: operators.f90:34
subroutine, public opgrad(ux, uy, uz, u, coef, es, ee)
Compute the weak gradient of a scalar field, i.e. the gradient multiplied by the mass matrix.
Definition: operators.f90:171
subroutine, public cdtp(dtx, x, dr, ds, dt, coef, es, ee)
Apply D^T to a scalar field, where D is the derivative matrix.
Definition: operators.f90:230
Krylov preconditioner.
Definition: precon.f90:34
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
Compound scheme for the advection and diffusion operators in a transport equation.
Base type for a matrix-vector product providing .
Definition: ax.f90:43
A list of boundary conditions.
Definition: bc.f90:104
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
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:66
Defines a canonical Krylov preconditioner.
Definition: precon.f90:40
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...