Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
adv_adjoint_dealias.f90
Go to the documentation of this file.
1
34!
36module adv_lin_dealias
37 use advection_adjoint, only: advection_adjoint_t
38 use num_types, only: rp
39 use math, only: vdot3, sub2
40 use space, only: space_t, gl
41 use field, only: field_t
42 use coefs, only: coef_t
43 use neko_config, only: neko_bcknd_device, neko_bcknd_sx, neko_bcknd_xsmm, &
44 neko_bcknd_opencl, neko_bcknd_cuda, neko_bcknd_hip
45 use operators, only: opgrad, cdtp
46 use interpolation, only: interpolator_t
47 use device_math, only: device_vdot3, device_sub2, device_col3, device_add4
48 use device, only: device_map
49 use utils, only: neko_error
50 use, intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr
51 implicit none
52 private
53
55 type, public, extends(advection_adjoint_t) :: adv_lin_dealias_t
57 type(coef_t) :: coef_gl
59 type(coef_t), pointer :: coef_gll
61 type(interpolator_t) :: gll_to_gl
63 type(space_t) :: xh_gl
65 type(space_t), pointer :: xh_gll
66 real(kind=rp), allocatable :: temp(:), tbf(:)
68 real(kind=rp), allocatable :: tx(:), ty(:), tz(:)
69 real(kind=rp), allocatable :: vr(:), vs(:), vt(:)
71 type(c_ptr) :: temp_d = c_null_ptr
73 type(c_ptr) :: tbf_d = c_null_ptr
75 type(c_ptr) :: tx_d = c_null_ptr
77 type(c_ptr) :: ty_d = c_null_ptr
79 type(c_ptr) :: tz_d = c_null_ptr
81 type(c_ptr) :: vr_d = c_null_ptr
83 type(c_ptr) :: vs_d = c_null_ptr
85 type(c_ptr) :: vt_d = c_null_ptr
86
87 real(kind=rp), allocatable :: txb(:), tyb(:), tzb(:)
89 type(c_ptr) :: txb_d = c_null_ptr
91 type(c_ptr) :: tyb_d = c_null_ptr
93 type(c_ptr) :: tzb_d = c_null_ptr
94
95 real(kind=rp), allocatable :: duxb(:), duyb(:), duzb(:)
96 real(kind=rp), allocatable :: dvxb(:), dvyb(:), dvzb(:)
97 real(kind=rp), allocatable :: dwxb(:), dwyb(:), dwzb(:)
99 type(c_ptr) :: duxb_d = c_null_ptr
101 type(c_ptr) :: duyb_d = c_null_ptr
103 type(c_ptr) :: duzb_d = c_null_ptr
105 type(c_ptr) :: dvxb_d = c_null_ptr
107 type(c_ptr) :: dvyb_d = c_null_ptr
109 type(c_ptr) :: dvzb_d = c_null_ptr
111 type(c_ptr) :: dwxb_d = c_null_ptr
113 type(c_ptr) :: dwyb_d = c_null_ptr
115 type(c_ptr) :: dwzb_d = c_null_ptr
116
117 contains
121 procedure, pass(this) :: compute_linear => &
122 compute_linear_advection_dealias
128 procedure, pass(this) :: compute_adjoint => &
129 compute_adjoint_advection_dealias
131 ! If one integrates by parts, this essentially switches sign and adds some
132 ! boundary terms.
133 ! We keep the differential operator on the test function
134 procedure, pass(this) :: compute_adjoint_scalar => &
135 compute_adjoint_scalar_advection_dealias
136 ! NOTE
137 ! This linearized advection term is the same as a normal advection term
138 ! so not sure what to do here...
140 procedure, pass(this) :: init => init_dealias
142 procedure, pass(this) :: free => free_dealias
143 end type adv_lin_dealias_t
144
145contains
146
151 subroutine init_dealias(this, lxd, coef)
152 class(adv_lin_dealias_t), target, intent(inout) :: this
153 integer, intent(in) :: lxd
154 type(coef_t), intent(inout), target :: coef
155 integer :: nel, n_GL, n
156
157 call this%Xh_GL%init(gl, lxd, lxd, lxd)
158 this%Xh_GLL => coef%Xh
159 this%coef_GLL => coef
160 call this%GLL_to_GL%init(this%Xh_GL, this%Xh_GLL)
161
162 call this%coef_GL%init(this%Xh_GL, coef%msh)
163
164 nel = coef%msh%nelv
165 n_gl = nel*this%Xh_GL%lxyz
166 n = nel*coef%Xh%lxyz
167 call this%GLL_to_GL%map(this%coef_GL%drdx, coef%drdx, nel, this%Xh_GL)
168 call this%GLL_to_GL%map(this%coef_GL%dsdx, coef%dsdx, nel, this%Xh_GL)
169 call this%GLL_to_GL%map(this%coef_GL%dtdx, coef%dtdx, nel, this%Xh_GL)
170 call this%GLL_to_GL%map(this%coef_GL%drdy, coef%drdy, nel, this%Xh_GL)
171 call this%GLL_to_GL%map(this%coef_GL%dsdy, coef%dsdy, nel, this%Xh_GL)
172 call this%GLL_to_GL%map(this%coef_GL%dtdy, coef%dtdy, nel, this%Xh_GL)
173 call this%GLL_to_GL%map(this%coef_GL%drdz, coef%drdz, nel, this%Xh_GL)
174 call this%GLL_to_GL%map(this%coef_GL%dsdz, coef%dsdz, nel, this%Xh_GL)
175 call this%GLL_to_GL%map(this%coef_GL%dtdz, coef%dtdz, nel, this%Xh_GL)
176
177 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
178 (neko_bcknd_opencl .eq. 1) .or. (neko_bcknd_sx .eq. 1) .or. &
179 (neko_bcknd_xsmm .eq. 1)) then
180 allocate(this%temp(n_gl))
181 allocate(this%tbf(n_gl))
182 allocate(this%tx(n_gl))
183 allocate(this%ty(n_gl))
184 allocate(this%tz(n_gl))
185 allocate(this%vr(n_gl))
186 allocate(this%vs(n_gl))
187 allocate(this%vt(n_gl))
188 allocate(this%duxb(n_gl))
189 allocate(this%duyb(n_gl))
190 allocate(this%duzb(n_gl))
191 allocate(this%dvxb(n_gl))
192 allocate(this%dvyb(n_gl))
193 allocate(this%dvzb(n_gl))
194 allocate(this%dwxb(n_gl))
195 allocate(this%dwyb(n_gl))
196 allocate(this%dwzb(n_gl))
197 allocate(this%txb(n_gl))
198 allocate(this%tyb(n_gl))
199 allocate(this%tzb(n_gl))
200 end if
201
202 if (neko_bcknd_device .eq. 1) then
203 call device_map(this%temp, this%temp_d, n_gl)
204 call device_map(this%tbf, this%tbf_d, n_gl)
205 call device_map(this%tx, this%tx_d, n_gl)
206 call device_map(this%ty, this%ty_d, n_gl)
207 call device_map(this%tz, this%tz_d, n_gl)
208 call device_map(this%vr, this%vr_d, n_gl)
209 call device_map(this%vs, this%vs_d, n_gl)
210 call device_map(this%vt, this%vt_d, n_gl)
211 call device_map(this%duxb, this%duxb_d, n_gl)
212 call device_map(this%duyb, this%duyb_d, n_gl)
213 call device_map(this%duzb, this%duzb_d, n_gl)
214 call device_map(this%dvxb, this%dvxb_d, n_gl)
215 call device_map(this%dvyb, this%dvyb_d, n_gl)
216 call device_map(this%dvzb, this%dvzb_d, n_gl)
217 call device_map(this%dwxb, this%dwxb_d, n_gl)
218 call device_map(this%dwyb, this%dwyb_d, n_gl)
219 call device_map(this%dwzb, this%dwzb_d, n_gl)
220 call device_map(this%txb, this%txb_d, n_gl)
221 call device_map(this%tyb, this%tyb_d, n_gl)
222 call device_map(this%tzb, this%tzb_d, n_gl)
223 end if
224
225 end subroutine init_dealias
226
228 subroutine free_dealias(this)
229 class(adv_lin_dealias_t), intent(inout) :: this
230 end subroutine free_dealias
231
232
250 subroutine compute_adjoint_advection_dealias(this, vx, vy, vz, vxb, vyb, &
251 vzb, fx, fy, fz, Xh, coef, n)
252 !! HARRY added vxb etc for baseflow
253 implicit none
254 class(adv_lin_dealias_t), intent(inout) :: this
255 type(space_t), intent(inout) :: Xh
256 type(coef_t), intent(inout) :: coef
257 type(field_t), intent(inout) :: vx, vy, vz
258 type(field_t), intent(inout) :: vxb, vyb, vzb
259 type(field_t), intent(inout) :: fx, fy, fz
260 integer, intent(in) :: n
261 real(kind=rp), dimension(this%Xh_GL%lxyz) :: tfx, tfy, tfz
262
263 ! u and U_b on dealiased mesh (one element)
264 real(kind=rp), dimension(this%Xh_GL%lxyz) :: tx, ty, tz
265 real(kind=rp), dimension(this%Xh_GL%lxyz) :: txb, tyb, tzb
266
267 ! gradients of U_b on dealiased mesh (one element)
268 real(kind=rp), dimension(this%Xh_GL%lxyz) :: duxb, dvxb, dwxb
269 real(kind=rp), dimension(this%Xh_GL%lxyz) :: duyb, dvyb, dwyb
270 real(kind=rp), dimension(this%Xh_GL%lxyz) :: duzb, dvzb, dwzb
271
272 real(kind=rp), dimension(this%Xh_GLL%lxyz) :: tempx, tempy, tempz
273 integer :: e, i, idx, nel, n_GL
274
275 nel = coef%msh%nelv
276 n_gl = nel * this%Xh_GL%lxyz
277 associate(c_gl => this%coef_GL)
278
279 if (neko_bcknd_device .eq. 1) then
280 ! Map baseflow to GL
281 call this%GLL_to_GL%map(this%txb, vxb%x, nel, this%Xh_GL)
282 call this%GLL_to_GL%map(this%tyb, vyb%x, nel, this%Xh_GL)
283 call this%GLL_to_GL%map(this%tzb, vzb%x, nel, this%Xh_GL)
284
285 ! Map adjoint velocity to GL
286 call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
287 call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
288 call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
289
290
291 ! u . grad U_b^T
292 !-----------------------------
293 ! take all the gradients
294 call opgrad(this%duxb, this%duyb, this%duzb, this%txb, c_gl)
295 call opgrad(this%dvxb, this%dvyb, this%dvzb, this%tyb, c_gl)
296 call opgrad(this%dwxb, this%dwyb, this%dwzb, this%tzb, c_gl)
297
298 ! traspose and multiply
299 call device_vdot3(this%vr_d, this%tx_d, this%ty_d, this%tz_d, &
300 this%duxb_d, this%dvxb_d, this%dwxb_d, n_gl)
301 call this%GLL_to_GL%map(this%temp, this%vr, nel, this%Xh_GLL)
302 call device_sub2(fx%x_d, this%temp_d, n)
303
304 call device_vdot3(this%vr_d, this%tx_d, this%ty_d, this%tz_d, &
305 this%duyb_d, this%dvyb_d, this%dwyb_d, n_gl)
306 call this%GLL_to_GL%map(this%temp, this%vr, nel, this%Xh_GLL)
307 call device_sub2(fy%x_d, this%temp_d, n)
308
309 call device_vdot3(this%vr_d, this%tx_d, this%ty_d, this%tz_d, &
310 this%duzb_d, this%dvzb_d, this%dwzb_d, n_gl)
311 call this%GLL_to_GL%map(this%temp, this%vr, nel, this%Xh_GLL)
312 call device_sub2(fz%x_d, this%temp_d, n)
313
314 ! \int \grad v . U_b ^ u with ^ an outer product
315
316 ! (x)
317 ! duxb,duyb,duzb are temporary arrays
318 call device_col3(this%duxb_d, this%tx_d, this%txb_d, n_gl)
319 call device_col3(this%duyb_d, this%tx_d, this%tyb_d, n_gl)
320 call device_col3(this%duzb_d, this%tx_d, this%tzb_d, n_gl)
321
322 ! D^T
323 ! vr,vs,vt are temporary arrays
324 call cdtp(this%vr, this%duxb, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl)
325 call cdtp(this%vs, this%duyb, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl)
326 call cdtp(this%vt, this%duzb, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl)
327
328 ! reuse duxb as a temp
329 call device_add4(this%duxb_d, this%vr_d, this%vs_d, this%vt_d, n_gl)
330 call this%GLL_to_GL%map(this%temp, this%duxb, nel, this%Xh_GLL)
331 call device_sub2(fx%x_d, this%temp_d, n)
332
333
334 ! (y)
335 ! duxb,duyb,duzb are temporary arrays
336 call device_col3(this%duxb_d, this%ty_d, this%txb_d, n_gl)
337 call device_col3(this%duyb_d, this%ty_d, this%tyb_d, n_gl)
338 call device_col3(this%duzb_d, this%ty_d, this%tzb_d, n_gl)
339
340 ! D^T
341 ! vr,vs,vt are temporary arrays
342 call cdtp(this%vr, this%duxb, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl)
343 call cdtp(this%vs, this%duyb, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl)
344 call cdtp(this%vt, this%duzb, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl)
345
346 ! reuse duxb as a temp
347 call device_add4(this%duxb_d, this%vr_d, this%vs_d, this%vt_d, n_gl)
348 call this%GLL_to_GL%map(this%temp, this%duxb, nel, this%Xh_GLL)
349 call device_sub2(fy%x_d, this%temp_d, n)
350
351 ! (z)
352 ! duxb,duyb,duzb are temporary arrays
353 call device_col3(this%duxb_d, this%tz_d, this%txb_d, n_gl)
354 call device_col3(this%duyb_d, this%tz_d, this%tyb_d, n_gl)
355 call device_col3(this%duzb_d, this%tz_d, this%tzb_d, n_gl)
356
357 ! D^T
358 ! vr,vs,vt are temporary arrays
359 call cdtp(this%vr, this%duxb, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl)
360 call cdtp(this%vs, this%duyb, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl)
361 call cdtp(this%vt, this%duzb, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl)
362
363 ! reuse duxb as a temp
364 call device_add4(this%duxb_d, this%vr_d, this%vs_d, this%vt_d, n_gl)
365 call this%GLL_to_GL%map(this%temp, this%duxb, nel, this%Xh_GLL)
366 call device_sub2(fz%x_d, this%temp_d, n)
367 else if ((neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
368 !TODO
369 else
370
371
372
373 do e = 1, coef%msh%nelv
374 ! Map baseflow to GL
375 call this%GLL_to_GL%map(txb, vxb%x(1,1,1,e), 1, this%Xh_GL)
376 call this%GLL_to_GL%map(tyb, vyb%x(1,1,1,e), 1, this%Xh_GL)
377 call this%GLL_to_GL%map(tzb, vzb%x(1,1,1,e), 1, this%Xh_GL)
378
379 ! Map adjoint velocity to GL
380 call this%GLL_to_GL%map(tx, vx%x(1,1,1,e), 1, this%Xh_GL)
381 call this%GLL_to_GL%map(ty, vy%x(1,1,1,e), 1, this%Xh_GL)
382 call this%GLL_to_GL%map(tz, vz%x(1,1,1,e), 1, this%Xh_GL)
383
384
385 ! u . grad U_b^T
386 !-----------------------------
387 call opgrad(duxb, duyb, duzb, txb, c_gl, e, e)
388 call opgrad(dvxb, dvyb, dvzb, tyb, c_gl, e, e)
389 call opgrad(dwxb, dwyb, dwzb, tzb, c_gl, e, e)
390
391 ! traspose and multiply
392 do i = 1, this%Xh_GL%lxyz
393 tfx(i) = tx(i)*duxb(i) + ty(i)*dvxb(i) + tz(i)*dwxb(i)
394 tfy(i) = tx(i)*duyb(i) + ty(i)*dvyb(i) + tz(i)*dwyb(i)
395 tfz(i) = tx(i)*duzb(i) + ty(i)*dvzb(i) + tz(i)*dwzb(i)
396 end do
397
398 ! map back to GLL
399 call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
400 call this%GLL_to_GL%map(tempy, tfy, 1, this%Xh_GLL)
401 call this%GLL_to_GL%map(tempz, tfz, 1, this%Xh_GLL)
402
403 ! accumulate
404 idx = (e-1)*this%Xh_GLL%lxyz+1
405 call sub2(fx%x(idx, 1, 1, 1), tempx, this%Xh_GLL%lxyz)
406 call sub2(fy%x(idx, 1, 1, 1), tempy, this%Xh_GLL%lxyz)
407 call sub2(fz%x(idx, 1, 1, 1), tempz, this%Xh_GLL%lxyz)
408
409 ! (x)
410 do i = 1, this%Xh_GL%lxyz
411 duxb(i) = tx(i)*txb(i)
412 duyb(i) = tx(i)*tyb(i)
413 duzb(i) = tx(i)*tzb(i)
414 end do
415
416 ! D^T
417 call cdtp(tfx, duxb, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl, e, e)
418 call cdtp(tfy, duyb, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl, e, e)
419 call cdtp(tfz, duzb, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl, e, e)
420
421 ! sum them
422 do i = 1, this%Xh_GL%lxyz
423 tfx(i) = tfx(i) + tfy(i) + tfz(i)
424 end do
425
426 ! map back to GLL
427 call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
428 call sub2(fx%x(idx, 1, 1, 1), tempx, this%Xh_GLL%lxyz)
429
430 ! (y)
431 do i = 1, this%Xh_GL%lxyz
432 duxb(i) = ty(i)*txb(i)
433 duyb(i) = ty(i)*tyb(i)
434 duzb(i) = ty(i)*tzb(i)
435 end do
436
437 ! D^T
438 call cdtp(tfx, duxb, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl, e, e)
439 call cdtp(tfy, duyb, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl, e, e)
440 call cdtp(tfz, duzb, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl, e, e)
441
442 ! sum them
443 do i = 1, this%Xh_GL%lxyz
444 tfx(i) = tfx(i) + tfy(i) + tfz(i)
445 end do
446
447 ! map back to GLL
448 call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
449 call sub2(fy%x(idx, 1, 1, 1), tempx, this%Xh_GLL%lxyz)
450
451 ! (z)
452 do i = 1, this%Xh_GL%lxyz
453 duxb(i) = tz(i)*txb(i)
454 duyb(i) = tz(i)*tyb(i)
455 duzb(i) = tz(i)*tzb(i)
456 end do
457 ! D^T
458 call cdtp(tfx, duxb, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl, e, e)
459 call cdtp(tfy, duyb, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl, e, e)
460 call cdtp(tfz, duzb, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl, e, e)
461
462 ! sum them
463 do i = 1, this%Xh_GL%lxyz
464 tfx(i) = tfx(i) + tfy(i) + tfz(i)
465 end do
466
467 ! map back to GLL
468 call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
469 call sub2(fz%x(idx, 1, 1, 1), tempx, this%Xh_GLL%lxyz)
470
471 end do
472
473
474 end if
475 end associate
476
477 end subroutine compute_adjoint_advection_dealias
478
495 subroutine compute_linear_advection_dealias(this, vx, vy, vz, vxb, vyb, vzb, &
496 fx, fy, fz, Xh, coef, n)
497 implicit none
498 class(adv_lin_dealias_t), intent(inout) :: this
499 type(space_t), intent(inout) :: Xh
500 type(coef_t), intent(inout) :: coef
501 type(field_t), intent(inout) :: vx, vy, vz
502 type(field_t), intent(inout) :: vxb, vyb, vzb
503 integer, intent(in) :: n
504 type(field_t), intent(inout) :: fx, fy, fz
505 real(kind=rp), dimension(this%Xh_GL%lxyz) :: tx, ty, tz
506 real(kind=rp), dimension(this%Xh_GL%lxyz) :: txb, tyb, tzb
507 real(kind=rp), dimension(this%Xh_GL%lxyz) :: tfx, tfy, tfz
508 real(kind=rp), dimension(this%Xh_GL%lxyz) :: vr, vs, vt
509 real(kind=rp), dimension(this%Xh_GLL%lxyz) :: tempx, tempy, tempz
510
511 integer :: e, i, idx, nel, n_GL
512 nel = coef%msh%nelv
513 n_gl = nel * this%Xh_GL%lxyz
514
515 !This is extremely primitive and unoptimized on the device //Karp
516 associate(c_gl => this%coef_GL)
517
518 if (neko_bcknd_device .eq. 1) then
519 ! Map baseflow to GL
520 call this%GLL_to_GL%map(this%txb, vxb%x, nel, this%Xh_GL)
521 call this%GLL_to_GL%map(this%tyb, vyb%x, nel, this%Xh_GL)
522 call this%GLL_to_GL%map(this%tzb, vzb%x, nel, this%Xh_GL)
523
524 ! Map perturbed velocity to GL
525 call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
526 call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
527 call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
528
529 ! u'.grad U
530 call opgrad(this%vr, this%vs, this%vt, this%txb, c_gl)
531 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
532 this%tx_d, this%ty_d, this%tz_d, n_gl)
533 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
534 call device_sub2(fx%x_d, this%temp_d, n)
535
536
537 call opgrad(this%vr, this%vs, this%vt, this%tyb, c_gl)
538 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
539 this%tx_d, this%ty_d, this%tz_d, n_gl)
540 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
541 call device_sub2(fy%x_d, this%temp_d, n)
542
543 call opgrad(this%vr, this%vs, this%vt, this%tzb, c_gl)
544 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
545 this%tx_d, this%ty_d, this%tz_d, n_gl)
546 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
547 call device_sub2(fz%x_d, this%temp_d, n)
548
549 ! U.grad u'
550 call opgrad(this%vr, this%vs, this%vt, this%tx, c_gl)
551 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
552 this%txb_d, this%tyb_d, this%tzb_d, n_gl)
553 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
554 call device_sub2(fx%x_d, this%temp_d, n)
555
556
557 call opgrad(this%vr, this%vs, this%vt, this%ty, c_gl)
558 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
559 this%txb_d, this%tyb_d, this%tzb_d, n_gl)
560 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
561 call device_sub2(fy%x_d, this%temp_d, n)
562
563 call opgrad(this%vr, this%vs, this%vt, this%tz, c_gl)
564 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
565 this%txb_d, this%tyb_d, this%tzb_d, n_gl)
566 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
567 call device_sub2(fz%x_d, this%temp_d, n)
568
569 else if ((neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
570 ! Map baseflow to GL
571 call this%GLL_to_GL%map(this%txb, vxb%x, nel, this%Xh_GL)
572 call this%GLL_to_GL%map(this%tyb, vyb%x, nel, this%Xh_GL)
573 call this%GLL_to_GL%map(this%tzb, vzb%x, nel, this%Xh_GL)
574
575 ! Map perturbed velocity to GL
576 call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
577 call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
578 call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
579
580 ! u'.grad U
581 call opgrad(this%vr, this%vs, this%vt, this%txb, c_gl)
582 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
583 this%tx, this%ty, this%tz, n_gl)
584 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
585 call sub2(fx%x, this%temp, n)
586
587
588 call opgrad(this%vr, this%vs, this%vt, this%tyb, c_gl)
589 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
590 this%tx, this%ty, this%tz, n_gl)
591 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
592 call sub2(fy%x, this%temp, n)
593
594 call opgrad(this%vr, this%vs, this%vt, this%tzb, c_gl)
595 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
596 this%tx, this%ty, this%tz, n_gl)
597 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
598 call sub2(fz%x, this%temp, n)
599
600 ! U.grad u'
601 call opgrad(this%vr, this%vs, this%vt, this%tx, c_gl)
602 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
603 this%txb, this%tyb, this%tzb, n_gl)
604 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
605 call sub2(fx%x, this%temp, n)
606
607
608 call opgrad(this%vr, this%vs, this%vt, this%ty, c_gl)
609 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
610 this%txb, this%tyb, this%tzb, n_gl)
611 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
612 call sub2(fy%x, this%temp, n)
613
614 call opgrad(this%vr, this%vs, this%vt, this%tz, c_gl)
615 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
616 this%txb, this%tyb, this%tzb, n_gl)
617 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
618 call sub2(fz%x, this%temp, n)
619 else
620
621 do e = 1, coef%msh%nelv
622 ! Map baseflow to GL
623 call this%GLL_to_GL%map(txb, vxb%x(1,1,1,e), 1, this%Xh_GL)
624 call this%GLL_to_GL%map(tyb, vyb%x(1,1,1,e), 1, this%Xh_GL)
625 call this%GLL_to_GL%map(tzb, vzb%x(1,1,1,e), 1, this%Xh_GL)
626 ! Map perturbed velocity to GL
627 call this%GLL_to_GL%map(tx, vx%x(1,1,1,e), 1, this%Xh_GL)
628 call this%GLL_to_GL%map(ty, vy%x(1,1,1,e), 1, this%Xh_GL)
629 call this%GLL_to_GL%map(tz, vz%x(1,1,1,e), 1, this%Xh_GL)
630
631 ! u'.grad U
632 call opgrad(vr, vs, vt, txb, c_gl, e, e)
633 do i = 1, this%Xh_GL%lxyz
634 tfx(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
635 end do
636
637 call opgrad(vr, vs, vt, tyb, c_gl, e, e)
638 do i = 1, this%Xh_GL%lxyz
639 tfy(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
640 end do
641
642 call opgrad(vr, vs, vt, tzb, c_gl, e, e)
643 do i = 1, this%Xh_GL%lxyz
644 tfz(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
645 end do
646
647 call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
648 call this%GLL_to_GL%map(tempy, tfy, 1, this%Xh_GLL)
649 call this%GLL_to_GL%map(tempz, tfz, 1, this%Xh_GLL)
650
651 idx = (e-1)*this%Xh_GLL%lxyz+1
652 call sub2(fx%x(idx, 1, 1, 1), tempx, this%Xh_GLL%lxyz)
653 call sub2(fy%x(idx, 1, 1, 1), tempy, this%Xh_GLL%lxyz)
654 call sub2(fz%x(idx, 1, 1, 1), tempz, this%Xh_GLL%lxyz)
655
656 ! U.grad u'
657 call opgrad(vr, vs, vt, tx, c_gl, e, e)
658 do i = 1, this%Xh_GL%lxyz
659 tfx(i) = txb(i)*vr(i) + tyb(i)*vs(i) + tzb(i)*vt(i)
660 end do
661
662 call opgrad(vr, vs, vt, ty, c_gl, e, e)
663 do i = 1, this%Xh_GL%lxyz
664 tfy(i) = txb(i)*vr(i) + tyb(i)*vs(i) + tzb(i)*vt(i)
665 end do
666
667 call opgrad(vr, vs, vt, tz, c_gl, e, e)
668 do i = 1, this%Xh_GL%lxyz
669 tfz(i) = txb(i)*vr(i) + tyb(i)*vs(i) + tzb(i)*vt(i)
670 end do
671
672 call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
673 call this%GLL_to_GL%map(tempy, tfy, 1, this%Xh_GLL)
674 call this%GLL_to_GL%map(tempz, tfz, 1, this%Xh_GLL)
675
676 idx = (e-1)*this%Xh_GLL%lxyz+1
677 call sub2(fx%x(idx, 1, 1, 1), tempx, this%Xh_GLL%lxyz)
678 call sub2(fy%x(idx, 1, 1, 1), tempy, this%Xh_GLL%lxyz)
679 call sub2(fz%x(idx, 1, 1, 1), tempz, this%Xh_GLL%lxyz)
680 end do
681 end if
682 end associate
683 end subroutine compute_linear_advection_dealias
684
699 subroutine compute_adjoint_scalar_advection_dealias(this, vxb, vyb, vzb, s, &
700 fs, Xh, coef, n, dt)
701 class(adv_lin_dealias_t), intent(inout) :: this
702 type(field_t), intent(inout) :: vxb, vyb, vzb
703 type(field_t), intent(inout) :: s
704 type(field_t), intent(inout) :: fs
705 type(space_t), intent(inout) :: Xh
706 type(coef_t), intent(inout) :: coef
707 integer, intent(in) :: n
708 real(kind=rp), intent(in), optional :: dt
709
710 real(kind=rp), dimension(this%Xh_GL%lxyz) :: vx_gl, vy_gl, vz_gl, s_gl
711 real(kind=rp), dimension(this%Xh_GL%lxyz) :: work1, work2, work3
712 real(kind=rp), dimension(this%Xh_GL%lxyz) :: w1, w2, w3
713 real(kind=rp), dimension(this%Xh_GL%lxyz) :: f_gl
714 integer :: e, i, idx, nel, n_GL
715 real(kind=rp), dimension(this%Xh_GLL%lxyz) :: temp
716
717 nel = coef%msh%nelv
718 n_gl = nel * this%Xh_GL%lxyz
719
720 associate(c_gl => this%coef_GL)
721 if (neko_bcknd_device .eq. 1) then
722 ! Map baseflow to GL
723 call this%GLL_to_GL%map(this%txb, vxb%x, nel, this%Xh_GL)
724 call this%GLL_to_GL%map(this%tyb, vyb%x, nel, this%Xh_GL)
725 call this%GLL_to_GL%map(this%tzb, vzb%x, nel, this%Xh_GL)
726
727 ! Map adjoint scalar to GL (use tx as adjoint scalar array)
728 call this%GLL_to_GL%map(this%tx, s%x, nel, this%Xh_GL)
729
730 ! Outer product (use duxb, duyb, duzb as temporary arrays)
731 call device_col3(this%duxb_d, this%tx_d, this%txb_d, n_gl)
732 call device_col3(this%duyb_d, this%tx_d, this%tyb_d, n_gl)
733 call device_col3(this%duzb_d, this%tx_d, this%tzb_d, n_gl)
734
735 ! D^T
736 ! vr,vs,vt are temporary arrays
737 call cdtp(this%vr, this%duxb, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl)
738 call cdtp(this%vs, this%duyb, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl)
739 call cdtp(this%vt, this%duzb, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl)
740
741 ! reuse duxb as a temp for summing them
742 call device_add4(this%duxb_d, this%vr_d, this%vs_d, this%vt_d, n_gl)
743
744 ! map back to GLL
745 call this%GLL_to_GL%map(this%temp, this%duxb, nel, this%Xh_GLL)
746
747 !apply
748 call device_sub2(fs%x_d, this%temp_d, n)
749
750
751 else if ((neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
752 call neko_error("Adjoint scalar not implemented for SX")
753 else
754 do e = 1, coef%msh%nelv
755 ! Map baseflow to GL
756 call this%GLL_to_GL%map(vx_gl, vxb%x(1,1,1,e), 1, this%Xh_GL)
757 call this%GLL_to_GL%map(vy_gl, vyb%x(1,1,1,e), 1, this%Xh_GL)
758 call this%GLL_to_GL%map(vz_gl, vzb%x(1,1,1,e), 1, this%Xh_GL)
759
760 ! Map passive scalar velocity to GL
761 call this%GLL_to_GL%map(s_gl, s%x(1,1,1,e), 1, this%Xh_GL)
762
763 do i = 1, this%Xh_GL%lxyz
764 work1(i) = s_gl(i)*vx_gl(i)
765 work2(i) = s_gl(i)*vy_gl(i)
766 work3(i) = s_gl(i)*vz_gl(i)
767 end do
768
769 ! D^T
770 call cdtp(w1, work1, c_gl%drdx, c_gl%dsdx, c_gl%dtdx, c_gl, e, e)
771 call cdtp(w2, work2, c_gl%drdy, c_gl%dsdy, c_gl%dtdy, c_gl, e, e)
772 call cdtp(w3, work3, c_gl%drdz, c_gl%dsdz, c_gl%dtdz, c_gl, e, e)
773
774 ! sum them
775 do i = 1, this%Xh_GL%lxyz
776 f_gl(i) = w1(i) + w2(i) + w3(i)
777 end do
778
779 ! map back to GLL
780 idx = (e-1)*this%Xh_GLL%lxyz+1
781 call this%GLL_to_GL%map(temp, f_gl, 1, this%Xh_GLL)
782 call sub2(fs%x(idx, 1, 1, 1), temp, this%Xh_GLL%lxyz)
783
784 end do
785
786 end if
787 end associate
788
789 end subroutine compute_adjoint_scalar_advection_dealias
790end module adv_lin_dealias
adv_lin_dealias
Subroutines to add advection terms to the RHS of a transport equation.
Definition adv_adjoint_dealias.f90:36
adv_lin_dealias::init_dealias
subroutine init_dealias(this, lxd, coef)
Constructor.
Definition adv_adjoint_dealias.f90:152
advection_adjoint
Subroutines to add advection terms to the RHS of a transport equation.
Definition advection_adjoint.f90:36
adv_lin_dealias::adv_lin_dealias_t
Type encapsulating advection routines with dealiasing.
Definition adv_adjoint_dealias.f90:55
advection_adjoint::advection_adjoint_t
Base abstract type for computing the advection operator.
Definition advection_adjoint.f90:48