Neko  0.9.0
A portable framework for high-order spectral element flow simulations
force_torque.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, dp, sp
38  use json_module, only : json_file
42  use field, only : field_t
43  use operators, only : curl
44  use case, only : case_t
47  use field_writer, only : field_writer_t
48  use coefs, only : coef_t
49  use operators, only : strain_rate
50  use vector, only : vector_t
51  use dirichlet, only : dirichlet_t
52  use drag_torque
53  use logger, only : log_size, neko_log
54  use comm
55  use math, only : masked_red_copy, cadd, glsum, vcross
58  use device
59 
60  implicit none
61  private
62 
65  type, public, extends(simulation_component_t) :: force_torque_t
67  type(field_t), pointer :: u
69  type(field_t), pointer :: v
71  type(field_t), pointer :: w
73  type(field_t), pointer :: p
74 
75  !Masked working arrays
76  type(vector_t) :: n1, n2, n3
77  type(vector_t) :: r1, r2, r3
78  type(vector_t) :: force1, force2, force3
79  type(vector_t) :: force4, force5, force6
80  type(vector_t) :: pmsk
81  type(vector_t) :: s11msk, s22msk, s33msk, s12msk, s13msk, s23msk
82  real(kind=rp) :: center(3) = 0.0_rp
83  real(kind=rp) :: scale
84  integer :: zone_id
85  character(len=20) :: zone_name
86  type(coef_t), pointer :: coef
87  type(dirichlet_t) :: bc
88  character(len=80) :: print_format
89 
90  contains
92  procedure, pass(this) :: init => force_torque_init_from_json
94  procedure, pass(this) :: init_from_attributes => &
97  procedure, pass(this) :: free => force_torque_free
99  procedure, pass(this) :: compute_ => force_torque_compute
100  end type force_torque_t
101 
102 contains
103 
105  subroutine force_torque_init_from_json(this, json, case)
106  class(force_torque_t), intent(inout) :: this
107  type(json_file), intent(inout) :: json
108  class(case_t), intent(inout), target :: case
109  integer :: zone_id
110  real(kind=rp), allocatable :: center(:)
111  character(len=:), allocatable :: zone_name
112  real(kind=rp) :: scale
113  logical :: long_print
114 
115  ! Add fields keyword to the json so that the field_writer picks it up.
116  ! Will also add fields to the registry.
117 
118  call this%init_base(json, case)
119 
120  call json_get(json, 'zone_id', zone_id)
121  call json_get_or_default(json, 'zone_name', zone_name, ' ')
122  call json_get_or_default(json, 'scale', scale, 1.0_rp)
123  call json_get_or_default(json, 'long_print', long_print, .false.)
124  call json_get(json, 'center', center)
125  call force_torque_init_from_attributes(this, zone_id, zone_name, &
126  center, scale, case%fluid%c_xh, &
127  long_print)
128  end subroutine force_torque_init_from_json
129 
131  subroutine force_torque_init_from_attributes(this, zone_id, zone_name, &
132  center, scale, coef, long_print)
133  class(force_torque_t), intent(inout) :: this
134  real(kind=rp), intent(in) :: center(3)
135  real(kind=rp), intent(in) :: scale
136  character(len=*), intent(in) :: zone_name
137  integer, intent(in) :: zone_id
138  type(coef_t), target, intent(in) :: coef
139  logical, intent(in) :: long_print
140  integer :: n_pts, glb_n_pts, ierr
141  character(len=1000) :: log_buf
142 
143  this%coef => coef
144  this%zone_id = zone_id
145  this%center = center
146  this%scale = scale
147  this%zone_name = zone_name
148  if (long_print ) then
149  this%print_format = '(I7,E20.10,E20.10,E20.10,E20.10,A)'
150  else
151  this%print_format = '(I7,E13.5,E13.5,E13.5,E13.5,A)'
152  end if
153 
154  this%u => neko_field_registry%get_field_by_name("u")
155  this%v => neko_field_registry%get_field_by_name("v")
156  this%w => neko_field_registry%get_field_by_name("w")
157  this%p => neko_field_registry%get_field_by_name("p")
158 
159 
160  call this%bc%init_base(this%coef)
161  call this%bc%mark_zone(this%case%msh%labeled_zones(this%zone_id))
162  call this%bc%finalize()
163  n_pts = this%bc%msk(0)
164  call this%n1%init(n_pts)
165  call this%n2%init(n_pts)
166  call this%n3%init(n_pts)
167  call this%r1%init(n_pts)
168  call this%r2%init(n_pts)
169  call this%r3%init(n_pts)
170  call this%force1%init(n_pts)
171  call this%force2%init(n_pts)
172  call this%force3%init(n_pts)
173  call this%force4%init(n_pts)
174  call this%force5%init(n_pts)
175  call this%force6%init(n_pts)
176  call this%s11msk%init(n_pts)
177  call this%s22msk%init(n_pts)
178  call this%s33msk%init(n_pts)
179  call this%s12msk%init(n_pts)
180  call this%s13msk%init(n_pts)
181  call this%s23msk%init(n_pts)
182  call this%pmsk%init(n_pts)
183  call setup_normals(this%coef, this%bc%msk, this%bc%facet, &
184  this%n1%x, this%n2%x, this%n3%x, n_pts)
185  call masked_red_copy(this%r1%x, this%coef%dof%x, this%bc%msk, &
186  this%u%size(), n_pts)
187  call masked_red_copy(this%r2%x, this%coef%dof%y, this%bc%msk, &
188  this%u%size(), n_pts)
189  call masked_red_copy(this%r3%x, this%coef%dof%z, this%bc%msk, &
190  this%u%size(), n_pts)
191 
192  call mpi_allreduce(n_pts, glb_n_pts, 1, &
193  mpi_integer, mpi_sum, neko_comm, ierr)
194  ! Print some information
195  call neko_log%section('Force/torque calculation')
196  write(log_buf, '(A,I4,A,A)') 'Zone ', zone_id, ' ', trim(zone_name)
197  call neko_log%message(log_buf)
198  write(log_buf, '(A,I6, I6)') 'Global number of GLL points in zone: ', glb_n_pts
199  call neko_log%message(log_buf)
200  write(log_buf, '(A,E15.7,E15.7,E15.7)') 'Average of zone''s coordinates: ', &
201  glsum(this%r1%x, n_pts)/glb_n_pts, &
202  glsum(this%r2%x, n_pts)/glb_n_pts, &
203  glsum(this%r3%x, n_pts)/glb_n_pts
204  call neko_log%message(log_buf)
205  write(log_buf, '(A,E15.7,E15.7,E15.7)') 'Center for torque calculation: ', center
206  call neko_log%message(log_buf)
207  write(log_buf, '(A,E15.7)') 'Scale: ', scale
208  call neko_log%message(log_buf)
209  call neko_log%end_section()
210 
211  call cadd(this%r1%x,-center(1), n_pts)
212  call cadd(this%r2%x,-center(2), n_pts)
213  call cadd(this%r3%x,-center(3), n_pts)
214  if (neko_bcknd_device .eq. 1 .and. n_pts .gt. 0) then
215  call device_memcpy(this%n1%x, this%n1%x_d, n_pts, host_to_device, &
216  .false.)
217  call device_memcpy(this%n2%x, this%n2%x_d, n_pts, host_to_device, &
218  .false.)
219  call device_memcpy(this%n3%x, this%n3%x_d, n_pts, host_to_device, &
220  .true.)
221  call device_memcpy(this%r1%x, this%r1%x_d, n_pts, host_to_device, &
222  .false.)
223  call device_memcpy(this%r2%x, this%r2%x_d, n_pts, host_to_device, &
224  .false.)
225  call device_memcpy(this%r3%x, this%r3%x_d, n_pts, host_to_device, &
226  .true.)
227  end if
228 
229  end subroutine force_torque_init_from_attributes
230 
232  subroutine force_torque_free(this)
233  class(force_torque_t), intent(inout) :: this
234  call this%free_base()
235 
236  nullify(this%u)
237  nullify(this%v)
238  nullify(this%w)
239  nullify(this%p)
240  nullify(this%coef)
241  end subroutine force_torque_free
242 
246  subroutine force_torque_compute(this, t, tstep)
247  class(force_torque_t), intent(inout) :: this
248  real(kind=rp), intent(in) :: t
249  integer, intent(in) :: tstep
250  real(kind=rp) :: dgtq(12) = 0.0_rp
251  integer :: n_pts, temp_indices(6)
252  type(field_t), pointer :: s11, s22, s33, s12, s13, s23
253  character(len=1000) :: log_buf
254 
255  n_pts = this%bc%msk(0)
256 
257  call neko_scratch_registry%request_field(s11, temp_indices(1))
258  call neko_scratch_registry%request_field(s12, temp_indices(2))
259  call neko_scratch_registry%request_field(s13, temp_indices(3))
260  call neko_scratch_registry%request_field(s22, temp_indices(4))
261  call neko_scratch_registry%request_field(s23, temp_indices(5))
262  call neko_scratch_registry%request_field(s33, temp_indices(6))
263 
264  call strain_rate(s11%x, s22%x, s33%x, s12%x, &
265  s13%x, s23%x, this%u, this%v, this%w, this%coef)
266 
267  ! On the CPU we can actually just use the original subroutines...
268  if (neko_bcknd_device .eq. 0) then
269  call masked_red_copy(this%s11msk%x, s11%x, this%bc%msk, &
270  this%u%size(), n_pts)
271  call masked_red_copy(this%s22msk%x, s22%x, this%bc%msk, &
272  this%u%size(), n_pts)
273  call masked_red_copy(this%s33msk%x, s33%x, this%bc%msk, &
274  this%u%size(), n_pts)
275  call masked_red_copy(this%s12msk%x, s12%x, this%bc%msk, &
276  this%u%size(), n_pts)
277  call masked_red_copy(this%s13msk%x, s13%x, this%bc%msk, &
278  this%u%size(), n_pts)
279  call masked_red_copy(this%s23msk%x, s23%x, this%bc%msk, &
280  this%u%size(), n_pts)
281  call masked_red_copy(this%pmsk%x, this%p%x, this%bc%msk, &
282  this%u%size(), n_pts)
283  call calc_force_array(this%force1%x, this%force2%x, this%force3%x, &
284  this%force4%x, this%force5%x, this%force6%x, &
285  this%s11msk%x, &
286  this%s22msk%x, &
287  this%s33msk%x, &
288  this%s12msk%x, &
289  this%s13msk%x, &
290  this%s23msk%x, &
291  this%pmsk%x, &
292  this%n1%x, &
293  this%n2%x, &
294  this%n3%x, &
295  this%case%fluid%mu, &
296  n_pts)
297  dgtq(1) = glsum(this%force1%x, n_pts)
298  dgtq(2) = glsum(this%force2%x, n_pts)
299  dgtq(3) = glsum(this%force3%x, n_pts)
300  dgtq(4) = glsum(this%force4%x, n_pts)
301  dgtq(5) = glsum(this%force5%x, n_pts)
302  dgtq(6) = glsum(this%force6%x, n_pts)
303  !Overwriting masked s11, s22, s33 as they are no longer needed
304  this%s11msk = 0.0_rp
305  this%s22msk = 0.0_rp
306  this%s33msk = 0.0_rp
307  call vcross(this%s11msk%x, this%s22msk%x, this%s33msk%x, &
308  this%r1%x, this%r2%x, this%r3%x, &
309  this%force1%x, this%force2%x, this%force3%x, n_pts)
310 
311  dgtq(7) = glsum(this%s11msk%x, n_pts)
312  dgtq(8) = glsum(this%s22msk%x, n_pts)
313  dgtq(9) = glsum(this%s33msk%x, n_pts)
314  this%s11msk = 0.0_rp
315  this%s22msk = 0.0_rp
316  this%s33msk = 0.0_rp
317  call vcross(this%s11msk%x, this%s22msk%x, this%s33msk%x, &
318  this%r1%x, this%r2%x, this%r3%x, &
319  this%force4%x, this%force5%x, this%force6%x, n_pts)
320  dgtq(10) = glsum(this%s11msk%x, n_pts)
321  dgtq(11) = glsum(this%s22msk%x, n_pts)
322  dgtq(12) = glsum(this%s33msk%x, n_pts)
323  else
324  if (n_pts .gt. 0) then
325  call device_masked_red_copy(this%s11msk%x_d, s11%x_d, &
326  this%bc%msk_d, this%u%size(), n_pts)
327  call device_masked_red_copy(this%s22msk%x_d, s22%x_d, &
328  this%bc%msk_d, this%u%size(), n_pts)
329  call device_masked_red_copy(this%s33msk%x_d, s33%x_d, &
330  this%bc%msk_d, this%u%size(), n_pts)
331  call device_masked_red_copy(this%s12msk%x_d, s12%x_d, &
332  this%bc%msk_d, this%u%size(), n_pts)
333  call device_masked_red_copy(this%s13msk%x_d, s13%x_d, &
334  this%bc%msk_d, this%u%size(), n_pts)
335  call device_masked_red_copy(this%s23msk%x_d, s23%x_d, &
336  this%bc%msk_d, this%u%size(), n_pts)
337  call device_masked_red_copy(this%pmsk%x_d, this%p%x_d, &
338  this%bc%msk_d, this%u%size(), n_pts)
339 
340  call device_calc_force_array(this%force1%x_d, this%force2%x_d, &
341  this%force3%x_d, &
342  this%force4%x_d, &
343  this%force5%x_d, &
344  this%force6%x_d, &
345  this%s11msk%x_d, &
346  this%s22msk%x_d, &
347  this%s33msk%x_d, &
348  this%s12msk%x_d, &
349  this%s13msk%x_d, &
350  this%s23msk%x_d, &
351  this%pmsk%x_d, &
352  this%n1%x_d, &
353  this%n2%x_d, &
354  this%n3%x_d, &
355  this%case%fluid%mu, &
356  n_pts)
357  !Overwriting masked s11, s22, s33 as they are no longer needed
358  call device_vcross(this%s11msk%x_d, this%s22msk%x_d, &
359  this%s33msk%x_d, &
360  this%r1%x_d, this%r2%x_d, this%r3%x_d, &
361  this%force1%x_d, this%force2%x_d, &
362  this%force3%x_d, n_pts)
363  call device_vcross(this%s12msk%x_d,this%s13msk%x_d,this%s23msk%x_d, &
364  this%r1%x_d, this%r2%x_d, this%r3%x_d, &
365  this%force4%x_d, this%force5%x_d, this%force6%x_d, n_pts)
366  end if
367  dgtq(1) = device_glsum(this%force1%x_d, n_pts)
368  dgtq(2) = device_glsum(this%force2%x_d, n_pts)
369  dgtq(3) = device_glsum(this%force3%x_d, n_pts)
370  dgtq(4) = device_glsum(this%force4%x_d, n_pts)
371  dgtq(5) = device_glsum(this%force5%x_d, n_pts)
372  dgtq(6) = device_glsum(this%force6%x_d, n_pts)
373  dgtq(7) = device_glsum(this%s11msk%x_d, n_pts)
374  dgtq(8) = device_glsum(this%s22msk%x_d, n_pts)
375  dgtq(9) = device_glsum(this%s33msk%x_d, n_pts)
376  dgtq(10) = device_glsum(this%s12msk%x_d, n_pts)
377  dgtq(11) = device_glsum(this%s13msk%x_d, n_pts)
378  dgtq(12) = device_glsum(this%s23msk%x_d, n_pts)
379  end if
380  dgtq = this%scale*dgtq
381  write(log_buf,'(A, I4, A, A)') 'Force and torque on zone ', &
382  this%zone_id,' ', this%zone_name
383  call neko_log%message(log_buf)
384  write(log_buf,'(A)') &
385  'Time step, time, total force/torque, pressure, viscous, direction'
386  call neko_log%message(log_buf)
387  write(log_buf, this%print_format) &
388  tstep,t,dgtq(1)+dgtq(4),dgtq(1),dgtq(4),', forcex'
389  call neko_log%message(log_buf)
390  write(log_buf, this%print_format) &
391  tstep,t,dgtq(2)+dgtq(5),dgtq(2),dgtq(5),', forcey'
392  call neko_log%message(log_buf)
393  write(log_buf, this%print_format) &
394  tstep,t,dgtq(3)+dgtq(6),dgtq(3),dgtq(6),', forcez'
395  call neko_log%message(log_buf)
396  write(log_buf, this%print_format) &
397  tstep,t,dgtq(7)+dgtq(10),dgtq(7),dgtq(10),', torquex'
398  call neko_log%message(log_buf)
399  write(log_buf, this%print_format) &
400  tstep,t,dgtq(8)+dgtq(11),dgtq(8),dgtq(11),', torquey'
401  call neko_log%message(log_buf)
402  write(log_buf, this%print_format) &
403  tstep,t,dgtq(9)+dgtq(12),dgtq(9),dgtq(12),', torquez'
404  call neko_log%message(log_buf)
405  call neko_scratch_registry%relinquish_field(temp_indices)
406 
407  end subroutine force_torque_compute
408 
409 end module force_torque
Copy data between host and device (or device and device)
Definition: device.F90:51
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Definition: json_utils.f90:54
Retrieves a parameter by name or throws an error.
Definition: json_utils.f90:45
Defines a boundary condition.
Definition: bc.f90:34
Defines a simulation case.
Definition: case.f90:34
Coefficients.
Definition: coef.f90:34
Definition: comm.F90:1
type(mpi_comm) neko_comm
MPI communicator.
Definition: comm.F90:16
subroutine, public device_cadd(a_d, c, n)
Add a scalar to vector .
subroutine, public device_masked_red_copy(a_d, b_d, mask_d, n, m)
subroutine, public device_vcross(u1_d, u2_d, u3_d, v1_d, v2_d, v3_d, w1_d, w2_d, w3_d, n)
Compute a cross product (3-d version) assuming vector components etc.
real(kind=rp) function, public device_glsum(a_d, n)
Sum a vector of length n.
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
integer, parameter, public host_to_device
Definition: device.F90:47
Defines a dirichlet boundary condition.
Definition: dirichlet.f90:34
subroutine, public setup_normals(coef, mask, facets, n1, n2, n3, n_pts)
subroutine, public device_calc_force_array(force1, force2, force3, force4, force5, force6, s11, s22, s33, s12, s13, s23, p, n1, n2, n3, v, n_pts)
Calculate drag and torque from array of points.
subroutine, public calc_force_array(force1, force2, force3, force4, force5, force6, s11, s22, s33, s12, s13, s23, p, n1, n2, n3, v, n_pts)
Calculate drag and torque from array of points.
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Implements the field_writer_t type.
Defines a field.
Definition: field.f90:34
Implements fld_file_output_t.
Implements the force_torque_t type.
subroutine force_torque_init_from_attributes(this, zone_id, zone_name, center, scale, coef, long_print)
Actual constructor.
subroutine force_torque_free(this)
Destructor.
subroutine force_torque_compute(this, t, tstep)
Compute the force_torque field.
subroutine force_torque_init_from_json(this, json, case)
Constructor from json.
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
Logging routines.
Definition: log.f90:34
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
integer, parameter, public log_size
Definition: log.f90:42
Definition: math.f90:60
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition: math.f90:323
subroutine, public masked_red_copy(a, b, mask, n, m)
Copy a masked vector to reduced contigous vector .
Definition: math.f90:280
real(kind=rp) function, public glsum(a, n)
Sum a vector of length n.
Definition: math.f90:360
subroutine, public vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
Compute a cross product assuming vector components etc.
Definition: math.f90:514
integer, parameter, public dp
Definition: num_types.f90:9
integer, parameter, public sp
Definition: num_types.f90:8
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Operators.
Definition: operators.f90:34
subroutine, public curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
Definition: operators.f90:362
subroutine, public strain_rate(s11, s22, s33, s12, s13, s23, u, v, w, coef)
Compute the strain rate tensor, i.e 0.5 * du_i/dx_j + du_j/dx_i.
Definition: operators.f90:430
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.
Simulation components are objects that encapsulate functionality that can be fit to a particular comp...
Defines a vector.
Definition: vector.f90:34
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Generic Dirichlet boundary condition on .
Definition: dirichlet.f90:44
A simulation component that writes a 3d field to a file.
A simple output saving a list of fields to a .fld file.
A simulation component that computes the force_torque field. Added to the field registry as omega_x,...
Base abstract class for simulation components.