Neko  0.9.0
A portable framework for high-order spectral element flow simulations
spalding.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 !
33 !
35 module spalding
36  use field, only: field_t
37  use num_types, only : rp
38  use json_module, only : json_file
39  use dofmap, only : dofmap_t
40  use coefs, only : coef_t
41  use neko_config, only : neko_bcknd_device
42  use wall_model, only : wall_model_t
45  use logger, only : neko_log, neko_log_debug
46  use utils, only : neko_error
47 
48  implicit none
49  private
50 
53  type, public, extends(wall_model_t) :: spalding_t
55  real(kind=rp) :: kappa = 0.41_rp
57  real(kind=rp) :: b = 5.2_rp
58  contains
60  procedure, pass(this) :: init => spalding_init
62  procedure, pass(this) :: init_from_components => &
65  procedure, pass(this) :: free => spalding_free
67  procedure, pass(this) :: compute => spalding_compute
69  procedure, private, pass(this) :: solve
70  end type spalding_t
71 
72 contains
80  subroutine spalding_init(this, coef, msk, facet, nu, h_index, json)
81  class(spalding_t), intent(inout) :: this
82  type(coef_t), intent(in) :: coef
83  integer, intent(in) :: msk(:)
84  integer, intent(in) :: facet(:)
85  real(kind=rp), intent(in) :: nu
86  integer, intent(in) :: h_index
87  type(json_file), intent(inout) :: json
88  real(kind=rp) :: kappa, b
89 
90  call json_get_or_default(json, "kappa", kappa, 0.41_rp)
91  call json_get_or_default(json, "B", b, 5.2_rp)
92 
93  call this%init_from_components(coef, msk, facet, nu, h_index, kappa, b)
94  end subroutine spalding_init
95 
104  subroutine spalding_init_from_components(this, coef, msk, facet, nu, h_index,&
105  kappa, B)
106  class(spalding_t), intent(inout) :: this
107  type(coef_t), intent(in) :: coef
108  integer, intent(in) :: msk(:)
109  integer, intent(in) :: facet(:)
110  integer, intent(in) :: h_index
111  real(kind=rp), intent(in) :: nu
112  real(kind=rp), intent(in) :: kappa
113  real(kind=rp), intent(in) :: b
114 
115  if (neko_bcknd_device .eq. 1) then
116  call neko_error("Spalding's law is only available on the CPU backend.")
117  end if
118 
119  call this%init_base(coef, msk, facet, nu, h_index)
120 
121  this%kappa = kappa
122  this%B = b
123  end subroutine spalding_init_from_components
124 
125 
127  subroutine spalding_free(this)
128  class(spalding_t), intent(inout) :: this
129 
130  call this%free_base()
131 
132  end subroutine spalding_free
133 
137  subroutine spalding_compute(this, t, tstep)
138  class(spalding_t), intent(inout) :: this
139  real(kind=rp), intent(in) :: t
140  integer, intent(in) :: tstep
141  type(field_t), pointer :: u
142  type(field_t), pointer :: v
143  type(field_t), pointer :: w
144  integer :: i
145  real(kind=rp) :: ui, vi, wi, magu, utau, normu, guess
146 
147  u => neko_field_registry%get_field("u")
148  v => neko_field_registry%get_field("v")
149  w => neko_field_registry%get_field("w")
150 
151  do i=1, this%n_nodes
152  ! Sample the velocity
153  ui = u%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), this%ind_e(i))
154  vi = v%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), this%ind_e(i))
155  wi = w%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), this%ind_e(i))
156 
157  ! Project on tangential direction
158  normu = ui * this%n_x%x(i) + vi * this%n_y%x(i) + wi * this%n_z%x(i)
159 
160  ui = ui - normu * this%n_x%x(i)
161  vi = vi - normu * this%n_y%x(i)
162  wi = wi - normu * this%n_z%x(i)
163 
164  magu = sqrt(ui**2 + vi**2 + wi**2)
165 
166  ! Get initial guess for Newton solver
167  if (tstep .eq. 1) then
168  guess = sqrt(magu * this%nu / this%h%x(i))
169  else
170  guess = this%tau_x(i)**2 + this%tau_y(i)**2 + this%tau_z(i)**2
171  guess = sqrt(sqrt(guess))
172  end if
173 
174  utau = this%solve(magu, this%h%x(i), guess)
175 
176  ! Distribute according to the velocity vector
177  this%tau_x(i) = -utau**2 * ui / magu
178  this%tau_y(i) = -utau**2 * vi / magu
179  this%tau_z(i) = -utau**2 * wi / magu
180  end do
181 
182  end subroutine spalding_compute
183 
188  function solve(this, u, y, guess) result(utau)
189  class(spalding_t), intent(inout) :: this
190  real(kind=rp), intent(in) :: u
191  real(kind=rp), intent(in) :: y
192  real(kind=rp), intent(in) :: guess
193  real(kind=rp) :: yp, up, kappa, b, utau
194  real(kind=rp) :: error, f, df, old
195  integer :: niter, k, maxiter
196 
197  utau = guess
198  kappa = this%kappa
199  b = this%B
200 
201  maxiter = 100
202 
203  do k=1, maxiter
204  up = u / utau
205  yp = y * utau / this%nu
206  niter = k
207  old = utau
208 
209  ! Evaluate function and its derivative
210  f = (up + exp(-kappa*b)* &
211  (exp(kappa*up) - 1.0_rp - kappa*up - 0.5_rp*(kappa*up)**2 - &
212  1.0_rp/6*(kappa*up)**3) - yp)
213 
214  df = (-y / this%nu - u/utau**2 - kappa*up/utau*exp(-kappa*b) * &
215  (exp(kappa*up) - 1 - kappa*up - 0.5*(kappa*up)**2))
216 
217  ! Update solution
218  utau = utau - f / df
219 
220  error = abs((old - utau)/old)
221 
222  if (error < 1e-3) then
223  exit
224  endif
225 
226  enddo
227 
228  if ((niter .eq. maxiter) .and. (neko_log%level_ .eq. neko_log_debug)) then
229  write(*,*) "Newton not converged", error, f, utau, old, guess
230  end if
231 end function solve
232 
233 
234 end module spalding
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Definition: json_utils.f90:54
Coefficients.
Definition: coef.f90:34
Defines a mapping of the degrees of freedom.
Definition: dofmap.f90:35
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
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
Logging routines.
Definition: log.f90:34
integer, parameter, public neko_log_debug
Debug log level.
Definition: log.f90:73
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
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
Implements spalding_t.
Definition: spalding.f90:35
real(kind=rp) function solve(this, u, y, guess)
Newton solver for the algebraic equation defined by the law.
Definition: spalding.f90:189
subroutine spalding_compute(this, t, tstep)
Compute the wall shear stress.
Definition: spalding.f90:138
subroutine spalding_init(this, coef, msk, facet, nu, h_index, json)
Constructor from JSON.
Definition: spalding.f90:81
subroutine spalding_init_from_components(this, coef, msk, facet, nu, h_index, kappa, B)
Constructor from components.
Definition: spalding.f90:106
subroutine spalding_free(this)
Destructor for the spalding_t (base) class.
Definition: spalding.f90:128
Utilities.
Definition: utils.f90:35
Implements wall_model_t.
Definition: wall_model.f90:35
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Wall model based on Spalding's law of the wall. Reference: http://dx.doi.org/10.1115/1....
Definition: spalding.f90:53
Base abstract type for wall-stress models for wall-modelled LES.
Definition: wall_model.f90:54