38 use json_module,
only : json_file
81 type(
vector_t) :: s11msk, s22msk, s33msk, s12msk, s13msk, s23msk
82 real(kind=
rp) :: center(3) = 0.0_rp
83 real(kind=
rp) :: scale
85 character(len=20) :: zone_name
88 character(len=80) :: print_format
94 procedure, pass(this) :: init_from_attributes => &
107 type(json_file),
intent(inout) :: json
108 class(
case_t),
intent(inout),
target :: case
110 real(kind=
rp),
allocatable :: center(:)
111 character(len=:),
allocatable :: zone_name
112 real(kind=
rp) :: scale
113 logical :: long_print
118 call this%init_base(json,
case)
120 call json_get(json,
'zone_id', zone_id)
124 call json_get(json,
'center', center)
126 center, scale,
case%fluid%c_xh, &
132 center, scale, coef, long_print)
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
144 this%zone_id = zone_id
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)'
151 this%print_format =
'(I7,E13.5,E13.5,E13.5,E13.5,A)'
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)
184 this%n1%x, this%n2%x, this%n3%x, n_pts)
186 this%u%size(), n_pts)
188 this%u%size(), n_pts)
190 this%u%size(), n_pts)
192 call mpi_allreduce(n_pts, glb_n_pts, 1, &
195 call neko_log%section(
'Force/torque calculation')
196 write(log_buf,
'(A,I4,A,A)')
'Zone ', zone_id,
' ', trim(zone_name)
198 write(log_buf,
'(A,I6, I6)')
'Global number of GLL points in zone: ', glb_n_pts
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
205 write(log_buf,
'(A,E15.7,E15.7,E15.7)')
'Center for torque calculation: ', center
207 write(log_buf,
'(A,E15.7)')
'Scale: ', scale
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
234 call this%free_base()
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
255 n_pts = this%bc%msk(0)
265 s13%x, s23%x, this%u, this%v, this%w, this%coef)
268 if (neko_bcknd_device .eq. 0)
then
270 this%u%size(), n_pts)
272 this%u%size(), n_pts)
274 this%u%size(), n_pts)
276 this%u%size(), n_pts)
278 this%u%size(), n_pts)
280 this%u%size(), n_pts)
282 this%u%size(), n_pts)
284 this%force4%x, this%force5%x, this%force6%x, &
295 this%case%fluid%mu, &
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)
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)
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)
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)
324 if (n_pts .gt. 0)
then
326 this%bc%msk_d, this%u%size(), n_pts)
328 this%bc%msk_d, this%u%size(), n_pts)
330 this%bc%msk_d, this%u%size(), n_pts)
332 this%bc%msk_d, this%u%size(), n_pts)
334 this%bc%msk_d, this%u%size(), n_pts)
336 this%bc%msk_d, this%u%size(), n_pts)
338 this%bc%msk_d, this%u%size(), n_pts)
355 this%case%fluid%mu, &
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)
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
384 write(log_buf,
'(A)') &
385 'Time step, time, total force/torque, pressure, viscous, direction'
387 write(log_buf, this%print_format) &
388 tstep,t,dgtq(1)+dgtq(4),dgtq(1),dgtq(4),
', forcex'
390 write(log_buf, this%print_format) &
391 tstep,t,dgtq(2)+dgtq(5),dgtq(2),dgtq(5),
', forcey'
393 write(log_buf, this%print_format) &
394 tstep,t,dgtq(3)+dgtq(6),dgtq(3),dgtq(6),
', forcez'
396 write(log_buf, this%print_format) &
397 tstep,t,dgtq(7)+dgtq(10),dgtq(7),dgtq(10),
', torquex'
399 write(log_buf, this%print_format) &
400 tstep,t,dgtq(8)+dgtq(11),dgtq(8),dgtq(11),
', torquey'
402 write(log_buf, this%print_format) &
403 tstep,t,dgtq(9)+dgtq(12),dgtq(9),dgtq(12),
', torquez'
Copy data between host and device (or device and device)
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 boundary condition.
Defines a simulation case.
type(mpi_comm) neko_comm
MPI communicator.
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.
integer, parameter, public host_to_device
Defines a dirichlet boundary condition.
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.
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.
type(log_t), public neko_log
Global log stream.
integer, parameter, public log_size
subroutine, public cadd(a, s, n)
Add a scalar to vector .
subroutine, public masked_red_copy(a, b, mask, n, m)
Copy a masked vector to reduced contigous vector .
real(kind=rp) function, public glsum(a, n)
Sum a vector of length n.
subroutine, public vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
Compute a cross product assuming vector components etc.
integer, parameter, public dp
integer, parameter, public sp
integer, parameter, public rp
Global precision used in computations.
subroutine, public curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
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.
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...
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Generic Dirichlet boundary condition on .
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.