Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
adjoint_minimum_dissipation_source_term.f90
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!
34!
35!
36! If the objective function $\int |\nabla u|^2$,
37! the corresponding adjoint forcing is $ \int \nabla v \cdot \nabla u $ in weak
38! form.
40 use num_types, only: rp
41 use field_list, only: field_list_t
42 use json_module, only: json_file
43 use json_utils, only: json_get, json_get_or_default
44 use source_term, only: source_term_t
45 use coefs, only: coef_t
46 use neko_config, only: neko_bcknd_device
47 use utils, only: neko_error
48 use field, only: field_t
49 use field_math, only: field_subcol3, field_add2, field_add2s2
50 use json_module, only: json_file
51 use time_state, only: time_state_t
53 use simcomp_executor, only: neko_simcomps
54 use user_source_term, only: user_source_term_t
55 use num_types, only: rp
56 use field, only: field_t
57 use field_registry, only: neko_field_registry
58 use math, only: rzero, copy, chsign
59 use device_math, only: device_copy, device_cmult
60 use neko_config, only: neko_bcknd_device
61 use operators, only: curl
62 use scratch_registry, only: neko_scratch_registry
64 use point_zone, only: point_zone_t
65 implicit none
66 private
67
69 ! $\int \nabla v \cdot \nabla u $
70 type, public, extends(source_term_t) :: &
73 type(field_t), pointer :: u => null()
75 type(field_t), pointer :: v => null()
77 type(field_t), pointer :: w => null()
79 real(kind=rp) :: obj_scale
81 class(point_zone_t), pointer :: mask => null()
83 logical :: if_mask
84
85 contains
87 procedure, pass(this) :: init => &
88 adjoint_minimum_dissipation_source_term_init_from_json
90 procedure, pass(this) :: init_from_components => &
91 adjoint_minimum_dissipation_source_term_init_from_components
93 procedure, pass(this) :: free => &
94 adjoint_minimum_dissipation_source_term_free
96 procedure, pass(this) :: compute_ => &
97 adjoint_minimum_dissipation_source_term_compute
99
100contains
107 subroutine adjoint_minimum_dissipation_source_term_init_from_json(this, &
108 json, fields, coef, variable_name)
109 class(adjoint_minimum_dissipation_source_term_t), intent(inout) :: this
110 type(json_file), intent(inout) :: json
111 type(field_list_t), intent(in), target :: fields
112 type(coef_t), intent(in), target :: coef
113 character(len=*), intent(in) :: variable_name
114
115 ! we shouldn't be initializing this from JSON
116 ! maybe throw an error?
117
118
119 end subroutine adjoint_minimum_dissipation_source_term_init_from_json
120
129 subroutine adjoint_minimum_dissipation_source_term_init_from_components(this,&
130 f_x, f_y, f_z, &
131 u, v, w, obj_scale, &
132 mask, if_mask, &
133 coef)
134 class(adjoint_minimum_dissipation_source_term_t), intent(inout) :: this
135 type(field_t), pointer, intent(in) :: f_x, f_y, f_z
136 type(field_list_t) :: fields
137 type(coef_t) :: coef
138 real(kind=rp) :: start_time
139 real(kind=rp) :: end_time
140 real(kind=rp) :: obj_scale
141 type(field_t), intent(in), target :: u, v, w
142 class(point_zone_t), intent(in), target :: mask
143 logical :: if_mask
144
145 ! I wish you didn't need a start time and end time...
146 ! but I'm just going to set a super big number...
147 start_time = 0.0_rp
148 end_time = 100000000.0_rp
149
150 call this%free()
151
152 ! this is copying the fluid source term init
153 ! We package the fields for the source term to operate on in a field list.
154 call fields%init(3)
155 call fields%assign(1, f_x)
156 call fields%assign(2, f_y)
157 call fields%assign(3, f_z)
158
159 call this%init_base(fields, coef, start_time, end_time)
160
161 ! point everything in the correct places
162 this%u => u
163 this%v => v
164 this%w => w
165
166 this%obj_scale = obj_scale
167
168 this%if_mask = if_mask
169 if (this%if_mask) then
170 this%mask => mask
171 end if
172
173 end subroutine adjoint_minimum_dissipation_source_term_init_from_components
174
176 subroutine adjoint_minimum_dissipation_source_term_free(this)
177 class(adjoint_minimum_dissipation_source_term_t), intent(inout) :: this
178
179 call this%free_base()
180 end subroutine adjoint_minimum_dissipation_source_term_free
181
185 subroutine adjoint_minimum_dissipation_source_term_compute(this, time)
186 class(adjoint_minimum_dissipation_source_term_t), intent(inout) :: this
187 type(time_state_t), intent(in) :: time
188 type(field_t), pointer :: u, v, w
189 type(field_t), pointer :: fu, fv, fw
190 !type(field_t), pointer :: dudx, dudy, dudz
191 !type(field_t), pointer :: dvdx, dvdy, dvdz
192 !type(field_t), pointer :: dwdx, dwdy, dwdz
193 type(field_t), pointer :: wo1, wo2, wo3, wo4, wo5, wo6
194 type(field_t), pointer :: t1 , t2
195 integer :: temp_indices(8)
196 integer n
197
198
199 fu => this%fields%get_by_index(1)
200 fv => this%fields%get_by_index(2)
201 fw => this%fields%get_by_index(3)
202
203 n = fu%size()
204
205 ! fuck I'm not sure about this... I need a pen and paper
206 ! also there should be a way to pre-process this forcing term...
207 ! instead of recalculating it every time
208 u => this%u
209 v => this%v
210 w => this%w
211 call neko_scratch_registry%request_field(wo1, temp_indices(1))
212 call neko_scratch_registry%request_field(wo2, temp_indices(2))
213 call neko_scratch_registry%request_field(wo3, temp_indices(3))
214 call neko_scratch_registry%request_field(wo4, temp_indices(4))
215 call neko_scratch_registry%request_field(wo5, temp_indices(5))
216 call neko_scratch_registry%request_field(wo6, temp_indices(6))
217 call neko_scratch_registry%request_field(t1 , temp_indices(7))
218 call neko_scratch_registry%request_field(t2 , temp_indices(8))
219
220 ! TODO
221 ! ok we're computing gradients at every timestep... which is stupid...
222 ! BUT
223 ! if this was unsteady we would have to do this.
224
225 ! TODO
226 ! this is cheating a little bit...
227 ! in strong form, \nabla u . \nabla v => - v . \nabla^2 u + bdry
228 !
229 ! we can do this properly later in weak form, ideally using ax_helm or so
230 !
231 ! for now we'll work in strong form and ignore the bdry
232 ! and suffer the double derivative :/
233 !
234 ! in fact, we'll go even quicker and use
235 ! \nabla ^2 u = grad (div (u)) - curl ( curl (u ))
236 ! and assume divergence free u
237
238 ! => - \nabla ^2 u = curl ( curl (u ))
239
240 call curl(wo1, wo2, wo3, u, v, w, t1, t2, this%coef)
241 call curl(wo4, wo5, wo6, wo1, wo2, wo3, t1, t2, this%coef)
242
243 ! mask
244 if (this%if_mask) then
245 call mask_exterior_const(wo4, this%mask, 0.0_rp)
246 call mask_exterior_const(wo5, this%mask, 0.0_rp)
247 call mask_exterior_const(wo6, this%mask, 0.0_rp)
248 end if
249
250 call field_add2s2(fu, wo4, this%obj_scale)
251 call field_add2s2(fv, wo5, this%obj_scale)
252 call field_add2s2(fw, wo6, this%obj_scale)
253
254 ! don't worry... we'll write this MUCH cleaner in the final version
255
256 call neko_scratch_registry%relinquish_field(temp_indices)
257
258 end subroutine adjoint_minimum_dissipation_source_term_compute
259
Implements the adjoint_minimum_dissipation_source_term_t type.
Some common Masking operations we may need.
Definition mask_ops.f90:34
Implements the steady_simcomp_t type.
The steady_simcomp_t type is a simulation component that terminates a simulation when the normed diff...