Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
mma.f90
Go to the documentation of this file.
1!===========================================================================!
2! Method of Moving Asymptotes !
3! This implementation is based on the following documents: !
4! 1-https://people.kth.se/~krille/mmagcmma.pdf !
5! 2-https://people.kth.se/~krille/originalmma.pdf !
6! 2-https://comsolyar.com/wp-content/uploads/2020/03/gcmma.pdf !
7! ------------------------------------------------------------------------- !
8! !
9! This module solves the following original optimization problem: !
10! !
11! Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 ) !
12! subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m !
13! xmin_j <= x_j <= xmax_j, j = 1,...,n !
14! z >= 0, y_i >= 0, i = 1,...,m !
15! !
16! by first creating the following convex approximation of the original !
17! problem: !
18! !
19! Minimize sum_{j = 1,...,n} (p0j / (upp_j-x_j) + q0j / (x_j-low_j)) + !
20! a0*z + sum_i = 1,...,m(c_i*y_i + 0.5*d_i*y_i^2) !
21! subject to sum_{j = 1,...,n} (pij / (upp_j-x_j) + qij / (x_j-low_j)) + !
22! a_i*z + y_i <= b_i, i = 1,...,m !
23! xmin_j <= alpha_j <= x_j <= beta_j <= xmax_j j = 1,...,n !
24! y_i>=0 i = 1,...,m !
25! z>=0 !
26! !
27! note that based on eq(3.5) there should be r0 in the approximated problem !
28! however since it is just a constant added to a minimization problem, it !
29! is ignored. !
30! A primal-dual algorithm is then employed to solve the aproximated problem !
31! using interior point method. !
32!===========================================================================!
33
34module mma
35 ! Inclusions from Neko
36 use num_types, only: rp
37 use json_module, only: json_file
38 use json_utils, only: json_get_or_default
39 use vector, only: vector_t
40 use matrix, only: matrix_t
41 use mpi_f08, only: mpi_allreduce, mpi_integer, mpi_sum, mpi_comm_world, &
42 mpi_min, mpi_max
43 use comm, only: pe_rank, neko_comm, mpi_real_precision
44 use, intrinsic :: iso_fortran_env, only: stderr => error_unit
45 use utils, only: neko_error
46 use device, only: device_memcpy, host_to_device
47
48 implicit none
49 private
50
51 type, public :: mma_t
52 private
53 integer :: n, m, max_iter
54 real(kind=rp) :: a0, f0val, asyinit, asyincr, asydecr, epsimin, &
55 residumax, residunorm
56 type(vector_t) :: xold1, xold2, low, upp, alpha, beta, a, c, d, xmax, xmin
57 logical :: is_initialized = .false.
58 logical :: is_updated = .false.
59 character(len=:), allocatable :: bcknd
60
61 ! Internal dummy variables for MMA
62 type(vector_t) :: p0j, q0j
63 type(matrix_t) :: pij, qij
64 type(vector_t) :: bi
65
66 !---nesessary for KKT check after updating df0dx, fval, dfdx --------
67 real(kind=rp) :: z, zeta
68 type(vector_t) :: y, lambda, s, mu
69 type(vector_t) :: xsi, eta
70 contains
72 generic, public :: init => init_from_json, init_from_components
73 procedure, public, pass(this) :: init_from_json => mma_init_from_json
74 procedure, public, pass(this) :: init_from_components => &
75 mma_init_from_components
76 procedure, public, pass(this) :: free => mma_free
77 procedure, public, pass(this) :: get_n => mma_get_n
78 procedure, public, pass(this) :: get_m => mma_get_m
79 procedure, public, pass(this) :: get_residumax => mma_get_residumax
80 procedure, public, pass(this) :: get_residunorm => mma_get_residunorm
81 procedure, public, pass(this) :: get_max_iter => mma_get_max_iter
82
83 procedure, public, pass(this) :: kkt => mma_kkt
84 procedure, public, pass(this) :: update => mma_update
85
86 end type mma_t
87
88 interface
89 ! ======================================================================= !
90 ! interface for cpu backend module subroutines
91
93 module subroutine mma_update_cpu(this, iter, x, df0dx, fval, dfdx)
94 class(mma_t), intent(inout) :: this
95 integer, intent(in) :: iter
96 real(kind=rp), dimension(this%n), intent(inout) :: x
97 type(vector_t) :: df0dx, fval
98 type(matrix_t) :: dfdx
99 end subroutine mma_update_cpu
100
102 module subroutine mma_kkt_cpu(this, x, df0dx, fval, dfdx)
103 class(mma_t), intent(inout) :: this
104 real(kind=rp), dimension(this%n), intent(in) :: x
105 type(vector_t), intent(in) :: df0dx, fval
106 type(matrix_t), intent(in) :: dfdx
107 end subroutine mma_kkt_cpu
108
109 ! ======================================================================= !
110 ! interface for device backend module subroutines
111
113 module subroutine mma_update_device(this, iter, x, df0dx, fval, dfdx)
114 class(mma_t), intent(inout) :: this
115 integer, intent(in) :: iter
116 real(kind=rp), dimension(this%n), intent(inout) :: x
117 type(vector_t) :: df0dx, fval
118 type(matrix_t) :: dfdx
119 end subroutine mma_update_device
120
122 module subroutine mma_kkt_device(this, x, df0dx, fval, dfdx)
123 class(mma_t), intent(inout) :: this
124 real(kind=rp), dimension(this%n), intent(in) :: x
125 type(vector_t), intent(in) :: df0dx, fval
126 type(matrix_t), intent(in) :: dfdx
127 end subroutine mma_kkt_device
128
129 end interface
130
131contains
132
134 subroutine mma_init_from_json(this, x, n, m, json, scale, auto_scale)
135 ! ----------------------------------------------------- !
136 ! Initializing the mma object and all the parameters !
137 ! required for MMA method. (a_i, c_i, d_i, ...) !
138 ! x: the design varaibles(DV), n: number of DV, !
139 ! m: number of constraints !
140 ! !
141 ! Note that residumax & residunorm of the KKT conditions!
142 ! are initialized with 10^5. This is done to avoid !
143 ! unnecessary extera computation of KKT norms for the !
144 ! initial design. !
145 ! ----------------------------------------------------- !
146 class(mma_t), intent(inout) :: this
147 integer, intent(in) :: n, m
148 real(kind=rp), intent(in), dimension(n) :: x
149
150 type(json_file), intent(inout) :: json
151
152 ! Read the scaling info for fval and dfdx from json
153 real(kind=rp), intent(out) :: scale
154 logical, intent(out) :: auto_scale
155 ! -------------------------------------------------------------------!
156 ! Internal parameters for MMA !
157 ! Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 ) !
158 ! subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m !
159 ! xmin_j <= x_j <= xmax_j, j = 1,...,n !
160 ! z >= 0, y_i >= 0, i = 1,...,m !
161 ! -------------------------------------------------------------------!
162 real(kind=rp), dimension(n) :: xmax, xmin
163 real(kind=rp), dimension(m) :: a, c, d
164 character(len=:), allocatable :: bcknd
165
166 ! For reading the values from json and then set the value for the arrays
167 real(kind=rp) :: a0 , xmax_const, xmin_const, a_const, c_const, d_const
168
169 integer :: max_iter, n_global, ierr
170 real(kind=rp) :: epsimin, asyinit, asyincr, asydecr
171
172 call mpi_allreduce(n, n_global, 1, mpi_integer, &
173 mpi_sum, mpi_comm_world, ierr)
174
175 ! ------------------------------------------------------------------------ !
176 ! Assign defaults if nothing is parsed
177 ! based on the Cpp Code by Niels
178 call json_get_or_default(json, 'mma.epsimin', epsimin, &
179 1.0e-9_rp * sqrt(real(m + n_global, rp)))
180 call json_get_or_default(json, 'mma.max_iter', max_iter, 100)
181
182 ! Following parameters are set based on eq.3.8:--------
183 call json_get_or_default(json, 'mma.asyinit', asyinit, 0.5_rp)
184 call json_get_or_default(json, 'mma.asyincr', asyincr, 1.2_rp)
185 call json_get_or_default(json, 'mma.asydecr', asydecr, 0.7_rp)
186
187 call json_get_or_default(json, 'mma.backend', bcknd, 'cpu')
188
189 call json_get_or_default(json, 'mma.xmin', xmin_const, 0.0_rp)
190 call json_get_or_default(json, 'mma.xmax', xmax_const, 1.0_rp)
191 call json_get_or_default(json, 'mma.a0', a0, 1.0_rp)
192 call json_get_or_default(json, 'mma.a', a_const, 0.0_rp)
193 call json_get_or_default(json, 'mma.c', c_const, 100.0_rp)
194 call json_get_or_default(json, 'mma.d', d_const, 0.0_rp)
195
196 call json_get_or_default(json, 'mma.scale', scale, 10.0_rp)
197 call json_get_or_default(json, 'mma.auto_scale', auto_scale, .true.)
198
199 call json_get_or_default(json, 'mma.epsimin', epsimin, 1.0e-9_rp)
200 ! Initialize the MMA object with the parsed parameters
201 a = a_const
202 c = c_const
203 d = d_const
204 xmin = xmin_const
205 xmax = xmax_const
206 ! initializing the mma concrete type (mma_cpu_t or mma_device_t)
207 if (pe_rank .eq. 0) then
208 print *, "Initializing MMA backend to >>> ", bcknd
209 end if
210
211 ! ------------------------------------------------------------------------ !
212 ! Initialize the MMA object with the parameters read from json
213 ! call this%init(x, n, m, a0, a, c, d, xmin, xmax, &
214 ! max_iter, epsimin, asyinit, asyincr, asydecr, bcknd)
215 call this%init(x, n, m, a0, a, c, d, xmin, xmax, &
216 max_iter, epsimin, asyinit, asyincr, asydecr, bcknd)
217
218 end subroutine mma_init_from_json
219
221 subroutine mma_free(this)
222 class(mma_t), intent(inout) :: this
223 ! Deallocate the internal vectors
224 call this%xold1%free()
225 call this%xold2%free()
226 call this%alpha%free()
227 call this%beta%free()
228 call this%a%free()
229 call this%c%free()
230 call this%d%free()
231 call this%low%free()
232 call this%upp%free()
233 call this%xmax%free()
234 call this%xmin%free()
235 call this%p0j%free()
236 call this%q0j%free()
237 call this%bi%free()
238 call this%y%free()
239 call this%lambda%free()
240 call this%s%free()
241 call this%mu%free()
242 call this%xsi%free()
243 call this%eta%free()
244
245 ! Deallocate the internal dummy matrices
246 call this%pij%free()
247 call this%qij%free()
248
249 this%is_initialized = .false.
250 this%is_updated = .false.
251 end subroutine mma_free
252
254 subroutine mma_init_from_components(this, x, n, m, a0, a, c, d, xmin, xmax, &
255 max_iter, epsimin, asyinit, asyincr, asydecr, bcknd)
256 ! ----------------------------------------------------- !
257 ! Initializing the mma object and all the parameters !
258 ! required for MMA method. (a_i, c_i, d_i, ...) !
259 ! x: the design varaibles(DV), n: number of DV, !
260 ! m: number of constraints !
261 ! !
262 ! Note that residumax & residunorm of the KKT conditions!
263 ! are initialized with 10^5. This is done to avoid !
264 ! unnecessary extera computation of KKT norms for the !
265 ! initial design. !
266 ! ----------------------------------------------------- !
267 class(mma_t), intent(inout) :: this
268 integer, intent(in) :: n, m
269 real(kind=rp), intent(in), dimension(n) :: x
270 ! -------------------------------------------------------------------!
271 ! Internal parameters for MMA !
272 ! Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 ) !
273 ! subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m !
274 ! xmin_j <= x_j <= xmax_j, j = 1,...,n !
275 ! z >= 0, y_i >= 0, i = 1,...,m !
276 ! -------------------------------------------------------------------!
277 real(kind=rp), intent(in), dimension(n) :: xmax, xmin
278 real(kind=rp), intent(in), dimension(m) :: a, c, d
279 real(kind=rp), intent(in) :: a0
280 integer, intent(in), optional :: max_iter
281 real(kind=rp), intent(in), optional :: epsimin, asyinit, asyincr, asydecr
282 character(len=:), intent(in), allocatable :: bcknd
283
284 call this%free()
285
286 this%n = n
287 this%m = m
288
289 call this%xold1%init(n)
290 call this%xold2%init(n)
291 this%xold1%x = x
292 this%xold2%x = x
293
294 call this%alpha%init(n)
295 call this%beta%init(n)
296
297 call this%a%init(m)
298 call this%c%init(m)
299 call this%d%init(m)
300 call this%low%init(n)
301 call this%upp%init(n)
302 call this%xmax%init(n)
303 call this%xmin%init(n)
304
305 !internal dummy variables for MMA
306 call this%p0j%init(n)
307 call this%q0j%init(n)
308 call this%pij%init(m, n)
309 call this%qij%init(m, n)
310 call this%bi%init(m)
311
312 !---nesessary for KKT check after updating df0dx, fval, dfdx --------
313 call this%y%init(m)
314 call this%lambda%init(m)
315 call this%s%init(m)
316 call this%mu%init(m)
317 call this%xsi%init(n)
318 call this%eta%init(n)
319
320 this%a0 = a0
321 this%a%x = a
322 this%c%x = c
323 this%d%x = d
324
325 !setting the bounds for the design variable based on the problem
326 this%xmax%x = xmax
327 this%xmin%x = xmin
328
329 this%low%x(:) = minval(x)
330 this%upp%x(:) = maxval(x)
331
332 !setting KKT norms to a large number for the initial design
333 this%residumax = huge(0.0_rp)
334 this%residunorm = huge(0.0_rp)
335
336 ! Select backend type
337 select case (bcknd)
338 case ("cpu")
339 if (pe_rank == 0) then
340 print *, "MMA initialized with CPU backend!"
341 end if
342 case ("cuda")
343 ! upload all init values to device pointers
344 call device_memcpy(this%xold1%x, this%xold1%x_d, this%n, &
345 host_to_device, sync = .false.)
346 call device_memcpy(this%xold1%x, this%xold2%x_d, this%n, &
347 host_to_device, sync = .false.)
348 call device_memcpy(this%a%x, this%a%x_d, this%m, host_to_device, &
349 sync = .false.)
350 call device_memcpy(this%c%x, this%c%x_d, this%m, host_to_device, &
351 sync = .false.)
352 call device_memcpy(this%d%x, this%d%x_d, this%m, host_to_device, &
353 sync = .false.)
354 call device_memcpy(this%xmax%x, this%xmax%x_d, this%n, host_to_device, &
355 sync = .false.)
356 call device_memcpy(this%xmin%x, this%xmin%x_d, this%n, host_to_device, &
357 sync = .false.)
358 if (pe_rank == 0) then
359 print *, "MMA initialized with CUDA backend!"
360 end if
361 case default
362 call neko_error('Unknown backend in mma_init_attributes')
363 end select
364
365 ! ------------------------------------------------------------------------ !
366 ! Assign defaults if nothing is parsed
367
368 ! based on the Cpp Code by Niels
369 if (.not. present(epsimin)) this%epsimin = 1.0e-9_rp * sqrt(real(m + n, rp))
370 if (.not. present(max_iter)) this%max_iter = 100
371
372 ! Following parameters are set based on eq.3.8:--------
373 if (.not. present(asyinit)) this%asyinit = 0.5_rp
374 if (.not. present(asyincr)) this%asyincr = 1.2_rp
375 if (.not. present(asydecr)) this%asydecr = 0.7_rp
376
377 ! Assign values from inputs when present
378 if (present(max_iter)) this%max_iter = max_iter
379 if (present(epsimin)) this%epsimin = epsimin
380 if (present(asyinit)) this%asyinit = asyinit
381 if (present(asyincr)) this%asyincr = asyincr
382 if (present(asydecr)) this%asydecr = asydecr
383 this%bcknd = bcknd
384
385 if (pe_rank .eq. 0) then
386 print *, "MMA is initialized with a0 = ", a0, ", a = ", a, ", c = ", c, &
387 ", d = ", d, "epsimin = ", this%epsimin
388 end if
389 !the object is correctly initialized
390 this%is_initialized = .true.
391 end subroutine mma_init_from_components
392
394 subroutine mma_update(this, iter, x, df0dx, fval, dfdx)
395 class(mma_t), intent(inout) :: this
396 integer, intent(in) :: iter
397 real(kind=rp), dimension(this%n), intent(inout) :: x
398 type(vector_t) :: df0dx, fval
399 type(matrix_t) :: dfdx
400
401 ! Select backend type
402 select case (this%bcknd )
403 case ("cpu")
404 call mma_update_cpu(this, iter, x, df0dx, fval, dfdx)
405 case ("cuda")
406 call mma_update_device(this, iter, x, df0dx, fval, dfdx)
407 case default
408 call mma_update_cpu(this, iter, x, df0dx, fval, dfdx)
409 end select
410 end subroutine mma_update
411
413 subroutine mma_kkt(this, x, df0dx, fval, dfdx)
414 class(mma_t), intent(inout) :: this
415 real(kind=rp), dimension(this%n), intent(in) :: x
416 type(vector_t), intent(in) :: df0dx, fval
417 type(matrix_t), intent(in) :: dfdx
418
419 ! Select backend type
420 select case (this%bcknd )
421 case ("cpu")
422 call mma_kkt_cpu(this, x, df0dx, fval, dfdx)
423 case ("cuda")
424 call mma_kkt_device(this, x, df0dx, fval, dfdx)
425 case default
426 call mma_kkt_cpu(this,x, df0dx, fval, dfdx)
427 end select
428 end subroutine mma_kkt
429
430 ! ========================================================================== !
431 ! Getters and setters
432
434 pure function mma_get_n(this) result(n)
435 class(mma_t), intent(in) :: this
436 integer :: n
437 n = this%n
438 end function mma_get_n
439
441 pure function mma_get_m(this) result(m)
442 class(mma_t), intent(in) :: this
443 integer :: m
444 m = this%m
445 end function mma_get_m
446
448 pure function mma_get_residumax(this) result(residumax)
449 class(mma_t), intent(in) :: this
450 real(kind=rp) :: residumax
451 residumax = this%residumax
452 end function mma_get_residumax
453
455 pure function mma_get_residunorm(this) result(residunorm)
456 class(mma_t), intent(in) :: this
457 real(kind=rp) :: residunorm
458 residunorm = this%residunorm
459 end function mma_get_residunorm
460
462 pure function mma_get_max_iter(this) result(max_iter_value)
463 class(mma_t), intent(in) :: this
464 integer :: max_iter_value
465 max_iter_value = this%max_iter
466 end function mma_get_max_iter
467
468end module mma
469
mma
Definition mma.f90:34
mma::mma_get_residumax
pure real(kind=rp) function mma_get_residumax(this)
Get L^{inf} norm (Max Norm) of the KKT conditions.
Definition mma.f90:449
mma::mma_init_from_json
subroutine mma_init_from_json(this, x, n, m, json, scale, auto_scale)
CPU update function, runs one iteration of MMA.
Definition mma.f90:135
mma::mma_kkt
subroutine mma_kkt(this, x, df0dx, fval, dfdx)
Call the KKT ckeck function based on the backend.
Definition mma.f90:414
mma::mma_init_from_components
subroutine mma_init_from_components(this, x, n, m, a0, a, c, d, xmin, xmax, max_iter, epsimin, asyinit, asyincr, asydecr, bcknd)
Initialize the mma object based on the attributes from the json file.
Definition mma.f90:256
mma::mma_get_m
pure integer function mma_get_m(this)
Get the number of constriants.
Definition mma.f90:442
mma::mma_get_max_iter
pure integer function mma_get_max_iter(this)
Get the maximum number of iterations for the mma_subsolve inner loop.
Definition mma.f90:463
mma::mma_get_residunorm
pure real(kind=rp) function mma_get_residunorm(this)
Get L^{2} norm (Euclidean Norm) of the KKT conditions.
Definition mma.f90:456
mma::mma_get_n
pure integer function mma_get_n(this)
Get the number of design variables (nloc)
Definition mma.f90:435
mma::mma_update
subroutine mma_update(this, iter, x, df0dx, fval, dfdx)
Call the update function based on the backend.
Definition mma.f90:395
mma::mma_free
subroutine mma_free(this)
Deallocate mma object.
Definition mma.f90:222
mma::mma_t
Definition mma.f90:51