OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
hm_read_visc_prony.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| hm_read_visc_prony ../starter/source/materials/visc/hm_read_visc_prony.F
25!||--- called by ------------------------------------------------------
26!|| hm_read_visc ../starter/source/materials/visc/hm_read_visc.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| hm_get_float_array_index ../starter/source/devtools/hm_reader/hm_get_float_array_index.F
30!|| hm_get_floatv ../starter/source/devtools/hm_reader/hm_get_floatv.F
31!|| hm_get_floatv_dim ../starter/source/devtools/hm_reader/hm_get_floatv_dim.F
32!|| hm_get_intv ../starter/source/devtools/hm_reader/hm_get_intv.F
33!|| hm_option_is_encrypted ../starter/source/devtools/hm_reader/hm_option_is_encrypted.f
34!|| lm_least_square_prony ../starter/source/materials/visc/hm_read_visc_prony.F
35!|| lm_least_square_prony_2 ../starter/source/materials/visc/hm_read_visc_prony.F
36!||--- uses -----------------------------------------------------
37!|| hm_option_read_mod ../starter/share/modules1/hm_option_read_mod.F
38!|| message_mod ../starter/share/message_module/message_mod.F
39!|| submodel_mod ../starter/share/modules1/submodel_mod.F
40!|| table_mod ../starter/share/modules1/table_mod.f
41!||====================================================================
42 SUBROUTINE hm_read_visc_prony(VISC ,IVISC ,NTABLE ,TABLE ,MAT_ID ,
43 . UNITAB ,LSUBMODEL)
44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE visc_param_mod
48 USE unitab_mod
49 USE message_mod
50 USE submodel_mod
51 USE hm_option_read_mod
52 USE table_mod
53C-----------------------------------------------
54C I m p l i c i t T y p e s
55C-----------------------------------------------
56#include "implicit_f.inc"
57C-----------------------------------------------
58C C o m m o n B l o c k s
59C-----------------------------------------------
60#include "units_c.inc"
61C-----------------------------------------------
62C D u m m y A r g u m e n t s
63C-----------------------------------------------
64 INTEGER ,INTENT(IN) :: IVISC
65 INTEGER ,INTENT(IN) :: NTABLE
66 INTEGER ,INTENT(IN) :: MAT_ID
67 TYPE (VISC_PARAM_) ,INTENT(INOUT) :: VISC
68 TYPE (UNIT_TYPE_) ,INTENT(IN) :: UNITAB
69 TYPE (SUBMODEL_DATA),INTENT(IN) :: LSUBMODEL(*)
70 TYPE (TTABLE) ,INTENT(INOUT) :: TABLE(NTABLE)
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 INTEGER :: I,NUPARAM,NIPARAM,NPRONY,NUVAR,IFLAG,IMOD,ITAB,ISHAPE,
75 . fctid_g,fctid_k,fctid_gs,fctid_ks,fctid_gl,fctid_kl
76 my_real :: g(100),beta(100),k(100),betak(100)
77 my_real :: kv,costfg,costfk,derivg,derivk,ginfini,kinfini,
78 . xgscale,xkscale,xgscale_unit,xkscale_unit,
79 . ygscale,ykscale,ygscale_unit,ykscale_unit,
80 . xgs_scale,ygs_scale,xgs_scale_unit,ygs_scale_unit,
81 . xgl_scale,ygl_scale,xgl_scale_unit,ygl_scale_unit,
82 . xks_scale,yks_scale,xks_scale_unit,yks_scale_unit,
83 . xkl_scale,ykl_scale,xkl_scale_unit,ykl_scale_unit
84C
85 LOGICAL :: IS_AVAILABLE,IS_ENCRYPTED
86C=======================================================================
87 is_encrypted = .false.
88 is_available = .false.
89
90 CALL hm_option_is_encrypted(is_encrypted)
91C======================================
92 visc%ILAW = ivisc
93C
94 ! Table initialization
95 g(1:100) = zero
96 beta(1:100) = zero
97 k(1:100) = zero
98 betak(1:100) = zero
99
100C
101 !IFLAG = 0 this flag was dedicated to keep old fomrulation for version older than
102 ! 1st Card - Flags and prony order
103 CALL hm_get_intv ('Model_Order' ,nprony ,is_available,lsubmodel)
104 CALL hm_get_floatv ('MAT_K' ,kv ,is_available,lsubmodel,unitab)
105 CALL hm_get_intv ('MAT_Itab' ,itab ,is_available,lsubmodel)
106 IF (itab > 2) itab = 0
107 CALL hm_get_intv ('MAT_Ishape' ,ishape ,is_available,lsubmodel)
108C
109 IF (nprony == 0) CALL ancmsg(msgid=2026,msgtype=msgerror,
110 . anmode=aninfo_blind_1,i1=mat_id)
111C
112 IF (ishape > 1) ishape = 0
113C
114 ! 2nd Card - Input depending on flags
115 ! =======================================================================================
116 ! Tabulated inputs
117 ! =======================================================================================
118 ! -> Itab = 1 : functions of time
119 IF (itab == 1) THEN
120 !----------------------------------------------------------------------------------
121 ! -> Shear modulus VS time
122 !----------------------------------------------------------------------------------
123 CALL hm_get_intv ('Fct_G' ,fctid_g,is_available,lsubmodel)
124 CALL hm_get_floatv ('XGscale',xgscale,is_available,lsubmodel,unitab)
125 IF (xgscale == zero) THEN
126 CALL hm_get_floatv_dim('XGscale',xgscale_unit,is_available, lsubmodel, unitab)
127 xgscale = one * xgscale_unit
128 ENDIF
129 CALL hm_get_floatv ('YGscale',ygscale ,is_available,lsubmodel,unitab)
130 IF (ygscale == zero) THEN
131 CALL hm_get_floatv_dim('YGscale' ,ygscale_unit,is_available, lsubmodel, unitab)
132 ygscale = one * ygscale_unit
133 ENDIF
134 !----------------------------------------------------------------------------------
135 ! -> Bulk modulus VS time
136 !----------------------------------------------------------------------------------
137 CALL hm_get_intv ('Fct_K' ,fctid_k,is_available,lsubmodel)
138 CALL hm_get_floatv ('XKscale',xkscale,is_available,lsubmodel,unitab)
139 IF (xkscale == zero) THEN
140 CALL hm_get_floatv_dim('XKscale',xkscale_unit,is_available, lsubmodel, unitab)
141 xkscale = one * xkscale_unit
142 ENDIF
143 CALL hm_get_floatv ('YKscale',ykscale ,is_available,lsubmodel,unitab)
144 IF (ykscale == zero) THEN
145 CALL hm_get_floatv_dim('YKscale' ,ykscale_unit,is_available, lsubmodel, unitab)
146 ykscale = one * ykscale_unit
147 ENDIF
148 !----------------------------------------------------------------------------------------
149 ! Fit Prony series parameters to approximate tabulated functions by Least Squares method
150 !----------------------------------------------------------------------------------------
151 IF ((fctid_g > 0).AND.(nprony > 0)) THEN
152 CALL lm_least_square_prony(mat_id ,nprony ,table ,fctid_g ,xgscale ,ygscale ,
153 . g ,beta ,costfg ,derivg ,ishape ,ginfini )
154 ENDIF
155 IF ((fctid_k > 0).AND.(kv == zero).AND.(nprony > 0)) THEN
156 CALL lm_least_square_prony(mat_id ,nprony ,table ,fctid_k ,xkscale ,ykscale ,
157 . k ,betak ,costfk ,derivk ,ishape ,kinfini )
158 ENDIF
159 ! -> Itab = 2 ! loss and storage moduli as functions of frequencies
160 ELSEIF (itab == 2) THEN
161 !----------------------------------------------------------------------------------
162 ! -> Shear storage modulus VS frequency
163 !----------------------------------------------------------------------------------
164 CALL hm_get_intv ('Fct_Gs' ,fctid_gs ,is_available,lsubmodel)
165 CALL hm_get_floatv ('XGs_scale' ,xgs_scale ,is_available,lsubmodel,unitab)
166 IF (xgs_scale == zero) THEN
167 CALL hm_get_floatv_dim('XGs_scale',xgs_scale_unit,is_available, lsubmodel, unitab)
168 xgs_scale = one * xgs_scale_unit
169 ENDIF
170 CALL hm_get_floatv ('YGs_scale' ,ygs_scale ,is_available,lsubmodel,unitab)
171 IF (ygs_scale == zero) THEN
172 CALL hm_get_floatv_dim('YGs_scale',ygs_scale_unit,is_available, lsubmodel, unitab)
173 ygs_scale = one * ygs_scale_unit
174 ENDIF
175 !----------------------------------------------------------------------------------
176 ! -> Shear loss modulus VS frequency
177 !----------------------------------------------------------------------------------
178 CALL hm_get_intv ('Fct_Gl' ,fctid_gl ,is_available,lsubmodel)
179 CALL hm_get_floatv ('XGl_scale' ,xgl_scale ,is_available,lsubmodel,unitab)
180 IF (xgl_scale == zero) THEN
181 CALL hm_get_floatv_dim('XGl_scale',xgl_scale_unit,is_available, lsubmodel, unitab)
182 xgl_scale = one * xgl_scale_unit
183 ENDIF
184 CALL hm_get_floatv ('YGl_scale' ,ygl_scale ,is_available,lsubmodel,unitab)
185 IF (ygl_scale == zero) THEN
186 CALL hm_get_floatv_dim('YGl_scale',ygl_scale_unit,is_available, lsubmodel, unitab)
187 ygl_scale = one * ygl_scale_unit
188 ENDIF
189 !----------------------------------------------------------------------------------
190 ! -> Bulk storage modulus VS frequency
191 !----------------------------------------------------------------------------------
192 CALL hm_get_intv ('Fct_Ks' ,fctid_ks ,is_available,lsubmodel)
193 CALL hm_get_floatv ('XKs_scale' ,xks_scale ,is_available,lsubmodel,unitab)
194 IF (xks_scale == zero) THEN
195 CALL hm_get_floatv_dim('XKs_scale',xks_scale_unit,is_available, lsubmodel, unitab)
196 xks_scale = one * xks_scale_unit
197 ENDIF
198 CALL hm_get_floatv ('YKs_scale' ,yks_scale ,is_available,lsubmodel,unitab)
199 IF (yks_scale == zero) THEN
200 CALL hm_get_floatv_dim('YKs_scale',yks_scale_unit,is_available, lsubmodel, unitab)
201 yks_scale = one * yks_scale_unit
202 ENDIF
203 !----------------------------------------------------------------------------------
204 ! -> Bulk loss modulus VS frequency
205 !----------------------------------------------------------------------------------
206 CALL hm_get_intv ('Fct_Kl' ,fctid_kl ,is_available,lsubmodel)
207 CALL hm_get_floatv ('XKl_scale' ,xkl_scale ,is_available,lsubmodel,unitab)
208 IF (xkl_scale == zero) THEN
209 CALL hm_get_floatv_dim('XKl_scale',xkl_scale_unit,is_available, lsubmodel, unitab)
210 xkl_scale = one * xkl_scale_unit
211 ENDIF
212 CALL hm_get_floatv ('YKl_scale' ,ykl_scale ,is_available,lsubmodel,unitab)
213 IF (ykl_scale == zero) THEN
214 CALL hm_get_floatv_dim('YKl_scale',ykl_scale_unit,is_available, lsubmodel, unitab)
215 ykl_scale = one * ykl_scale_unit
216 ENDIF
217 !----------------------------------------------------------------------------------------
218 ! Fit Prony series parameters to approximate tabulated functions by Least Squares method
219 !----------------------------------------------------------------------------------------
220 IF ((fctid_gs > 0).AND.(fctid_gl > 0).AND.(nprony > 0)) THEN
221 CALL lm_least_square_prony_2(mat_id ,nprony ,table ,fctid_gs ,xgs_scale,ygs_scale,
222 . fctid_gl ,xgl_scale,ygl_scale,g ,beta ,costfg ,
223 . derivg ,ishape ,ginfini )
224 ENDIF
225 IF ((fctid_ks > 0).AND.(fctid_kl > 0).AND.(nprony > 0).AND.(kv == zero)) THEN
226 CALL lm_least_square_prony_2(mat_id ,nprony ,table ,fctid_ks ,xks_scale,yks_scale,
227 . fctid_kl ,xkl_scale,ykl_scale,k ,betak ,costfk ,
228 . derivk ,ishape ,kinfini )
229 ENDIF
230 ! =======================================================================================
231 ! Classical input
232 ! =======================================================================================
233 ! -> Itab = 0 ! classical input of prony series
234 ELSE
235 ishape = 0
236 IF(nprony > 0) THEN
237 DO i=1,nprony
238 CALL hm_get_float_array_index('Fport1' ,g(i) ,i,is_available,lsubmodel,unitab)
239 CALL hm_get_float_array_index('Fporp1' ,beta(i) ,i,is_available,lsubmodel,unitab)
240 CALL hm_get_float_array_index('Ki' ,k(i) ,i,is_available,lsubmodel,unitab)
241 CALL hm_get_float_array_index('Beta_ki',betak(i),i,is_available,lsubmodel,unitab)
242 ENDDO
243 ENDIF
244 ENDIF
245C
246 ! Taking into account the infinite value shape (+ 1 prony function with beta = 0)
247 IF ((itab /= 0) .AND. (ishape == 1)) THEN
248 nprony = nprony + 1
249 ENDIF
250c-------------------------------------------------
251c Storing parameters in UPARAM/IPARAM tables
252c-------------------------------------------------
253
254 nuvar = 7*nprony
255 niparam = 2
256 nuparam = 4*nprony + 1
257 ALLOCATE (visc%UPARAM(nuparam))
258 ALLOCATE (visc%IPARAM(niparam))
259 visc%NUVAR = nuvar
260 visc%NUPARAM = nuparam
261 visc%NIPARAM = niparam
262
263 imod = 0
264 visc%UPARAM(1) = kv
265 DO i=1,nprony
266 visc%UPARAM(1 + i) = g(i)
267 visc%UPARAM(1 + nprony + i) = beta(i)
268 visc%UPARAM(1 + 2*nprony + i) = k(i)
269 visc%UPARAM(1 + 3*nprony + i) = betak(i)
270 IF (k(i) > zero) imod = 1
271 ENDDO
272 visc%IPARAM(1) = nprony
273 visc%IPARAM(2) = imod
274c-----------------------------------------------------------------------
275c Output
276c-----------------------------------------------------------------------
277 IF (is_encrypted)THEN
278 WRITE(iout,'(5X,A,//)')'confidential data'
279 ELSE
280C
281 IF(NPRONY > 0) THEN
282 WRITE(IOUT,1000)
283 IF(IMOD == 0) THEN
284 WRITE(IOUT,1100) KV,NPRONY-ISHAPE
285 IF (ITAB > 0) THEN
286 WRITE(IOUT,1500) ITAB,ISHAPE
287 IF (ISHAPE == 1) WRITE(IOUT,3000) GINFINI
288 IF (ITAB == 1) THEN
289 WRITE(IOUT,1600) FctID_G,XGscale,YGscale,COSTFG,DERIVG
290 ELSEIF (ITAB == 2) THEN
291 WRITE(IOUT,2000) FctID_Gs,XGs_scale,YGs_scale,
292 . FctID_Gl,XGl_scale,YGl_scale,
293 . COSTFG,DERIVG
294 ENDIF
295 ENDIF
296 DO I=1,NPRONY-ISHAPE
297 WRITE(IOUT,1150) I
298 WRITE(IOUT,1200) G(I+ISHAPE),BETA(I+ISHAPE)
299 ENDDO
300 ELSE
301 WRITE(IOUT,1300) NPRONY-ISHAPE
302 IF (ITAB > 0) THEN
303 WRITE(IOUT,1500) ITAB,ISHAPE
304 IF (ITAB == 1) THEN
305 IF (ISHAPE == 1) WRITE(IOUT,3000) GINFINI
306 WRITE(IOUT,1600) FctID_G,XGscale,YGscale,COSTFG,DERIVG
307 IF (ISHAPE == 1) WRITE(IOUT,3500) KINFINI
308 WRITE(IOUT,1800) FctID_K,XKscale,YKscale,COSTFK,DERIVK
309 ELSEIF (ITAB == 2) THEN
310 IF (ISHAPE == 1) WRITE(IOUT,3000) GINFINI
311 WRITE(IOUT,2000) FctID_Gs,XGs_scale,YGs_scale,
312 . FctID_Gl,XGl_scale,YGl_scale,
313 . COSTFG,DERIVG
314 IF (ISHAPE == 1) WRITE(IOUT,3500) KINFINI
315 WRITE(IOUT,2500) FctID_Ks,XKs_scale,YKs_scale,
316 . FctID_Kl,XKl_scale,YKl_scale,
317 . COSTFK,DERIVK
318 ENDIF
319 ENDIF
320 DO I=1,NPRONY-ISHAPE
321 WRITE(IOUT,1150) I
322 WRITE(IOUT,1200) G(I+ISHAPE),BETA(I+ISHAPE)
323 WRITE(IOUT,1400) K(I+ISHAPE),BETAK(I+ISHAPE)
324 ENDDO
325 ENDIF
326 ENDIF
327 ENDIF
328C-----------
329 1000 FORMAT(
330 & 5X,' prony series model :' ,/,
331 & 5X,' --------------------- ' ,/)
332 1100 FORMAT(
333 & 5X,'bulk modulus for visco elastic material . . . . . . . . . . . . . . . =',1PG20.13/
334 & 5X,'order of prony series . . . . . . . . . . . . . . . . . . . . . . . . =',I10/)
335 1150 FORMAT(
336 & 5X,' ----------------------------------------------------------------------'/
337 & 5X,' parameters for prony FUNCTION number #',I10/
338 & 5X,' ----------------------------------------------------------------------'/)
339 1200 FORMAT(
340 & 5X,'shear relaxation g modulus . . . . . . . . . . . . . . . . . . . . . = '1PG20.13/
341 & 5X,'beta decay shear modulus . . . . . . . . . . . . . . . . . . . . . . =',1PG20.13)
342 1300 FORMAT(
343 & 5X,'order of prony series . . . . . . . . . . . . . . . . . . . . . . . . =',I10//)
344 1400 FORMAT(
345 & 5X,'bulk relaxation k modulus . . . . . . . . . . . . . . . . . . . . . . = '1PG20.13/
346 & 5X,'betak decay bulk modulus . . . . . . . . . . . . . . . . . . . . . . =',1PG20.13//)
347 1500 FORMAT(
348 & 5X,'tabulated prony series flag . . . . . . . . . . . . . . . . . . . . .=',I10/
349 & 5X,' itab=1 fitting from modulus vs time curves'/
350 & 5X,' itab=2 fitting from storage and loss moduli vs frequency curves'/
351 & 5X,'shape prony series flag . . . . . . . . . . . . . . . . . . . . . . .=',I10/
352 & 5X,' ishape=0 without infinite modulus (DEFAULT)'/
353 & 5X,' ishape=1 with infinite modulus'/)
354 1600 FORMAT(
355 & 5X,'least square fitting from shear modulus g function id. . . . . . . . .= 'I10/
356 & 5X,'time scale factor for shear modulus . . . . . . . . . . . . . . . . .= '1PG20.13/
357 & 5X,'scale factor for shear modulus . . . . . . . . . . . . . . . . . . . .= '1PG20.13/
358 & 5X,'final cost function value . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/
359 & 5X,'final derivative value . . . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/)
360 1800 FORMAT(
361 & 5X,'least square fitting from bulk modulus k function id . . . . . . . . .= 'I10/
362 & 5X,'time scale factor for bulk modulus . . . . . . . . . . . . . . . . . .= '1PG20.13/
363 & 5X,'scale factor for bulk modulus . . . . . . . . . . . . . . . . . . . .= '1PG20.13/
364 & 5X,'final cost function value . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/
365 & 5X,'final derivative value . . . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/)
366 2000 FORMAT(
367 & 5X,'least square fitting from storage shear modulus gl function id . . . .= 'I10/
368 & 5X,'frequency scale factor for storage shear modulus . . . . . . . . . . .= '1PG20.13/
369 & 5X,'scale factor for storage shear modulus . . . . . . . . . . . . . . . .= '1PG20.13/
370 & 5X,'least square fitting from loss shear modulus gs function id . . . . .= 'I10/
371 & 5X,'frequency scale factor for loss shear modulus . . . . . . . . . . . .= '1PG20.13/
372 & 5X,'scale factor for loss shear modulus function . . . . . . . . . . . . .= '1PG20.13/
373 & 5X,'final cost function value . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/
374 & 5X,'final derivative value . . . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/)
375 2500 FORMAT(
376 & 5X,'least square fitting from storage bulk modulus kl function id . . . .= 'I10/
377 & 5X,'frequency scale factor for storage bulk modulus . . . . . . . . . . .= '1PG20.13/
378 & 5X,'scale factor for storage bulk modulus . . . . . . . . . . . . . . . .= '1PG20.13/
379 & 5X,'least square fitting from loss bulk modulus gs function id . . . . . .= 'I10/
380 & 5X,'frequency scale factor for loss bulk modulus . . . . . . . . . . . . .= '1PG20.13/
381 & 5X,'scale factor for loss bulk modulus function . . . . . . . . . . . . .= '1PG20.13/
382 & 5X,'final cost function value . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/
383 & 5X,'final derivative value . . . . . . . . . . . . . . . . . . . . . . . .= '1PG20.13/)
384 3000 FORMAT(
385 & 5X,'shear modulus infinite value ginf . . . . . . . . . . . . . . . . . .= '1PG20.13/)
386 3500 FORMAT(
387 & 5X,'bulk modulus infinite value kinf . . . . . . . . . . . . . . . . . . .= '1PG20.13/)
388
389 RETURN
390 END
391
392!||====================================================================
393!|| lm_least_square_prony ../starter/source/materials/visc/hm_read_visc_prony.F
394!||--- called by ------------------------------------------------------
395!|| hm_read_visc_prony ../starter/source/materials/visc/hm_read_visc_prony.F
396!||--- calls -----------------------------------------------------
397!|| ancmsg ../starter/source/output/message/message.F
398!||--- uses -----------------------------------------------------
399!|| message_mod ../starter/share/message_module/message_mod.F
400!|| table_mod ../starter/share/modules1/table_mod.F
401!||====================================================================
402 SUBROUTINE LM_LEAST_SQUARE_PRONY(MAT_ID ,NPRONY ,TABLE ,FCT_ID ,XGSCALE ,YGSCALE ,
403 . G ,BETA ,COST ,DERIV ,ISHAPE ,GINFINI )
404C-----------------------------------------------
405C M o d u l e s
406C-----------------------------------------------
407 USE TABLE_MOD
408 USE MESSAGE_MOD
409C-----------------------------------------------
410C I m p l i c i t T y p e s
411C-----------------------------------------------
412#include "implicit_f.inc"
413#include "com04_c.inc"
414C-----------------------------------------------
415C Dummy arguments
416C-----------------------------------------------
417 INTEGER, INTENT(IN) :: NPRONY,FCT_ID,MAT_ID,ISHAPE
418 my_real
419 . XGSCALE,YGSCALE,COST,DERIV,GINFINI
420 my_real, DIMENSION(100) :: G,BETA
421 TYPE(TTABLE) TABLE(NTABLE)
422C-----------------------------------------------
423C Local arguments
424C-----------------------------------------------
425 INTEGER I,II,INFO,SIZE_TIME,NITER
426 INTEGER, DIMENSION(:), ALLOCATABLE :: IPIV
427 DOUBLE PRECISION ::
428 . RI2,RI2_OLD,MXVALT,MXVALY,LAMBDA,
429 . JACTDY_N2,RHO,NU,DENOM,XI_N2
430 DOUBLE PRECISION, PARAMETER :: EPS1 = 1.0D-8
431 DOUBLE PRECISION, PARAMETER :: EPS2 = 1.0D-10
432 DOUBLE PRECISION, PARAMETER :: ALPHA = 0.5D0
433 DOUBLE PRECISION, PARAMETER :: MXITER = 10000
434 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: XI,FXI,DELTAY,JACTDY,YI,TIME,
435 . XI_TR,FXI_TR,VECT,MU
436 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: JAC,JACTJAC,JACT,ID
437 LOGICAL :: CONV,FOUND
438C
439 ! Find the function and recover the experimental data
440 FOUND = .FALSE.
441 SIZE_TIME = 0
442 DO I = 1,NTABLE
443 IF (TABLE(I)%NOTABLE == FCT_ID) THEN
444 ! Number of data points
445 SIZE_TIME = SIZE(TABLE(I)%X(1)%VALUES)
446 ALLOCATE(YI(SIZE_TIME),TIME(SIZE_TIME))
447 ! Recovering data (time and modulus)
448 DO II = 1,SIZE_TIME
449 TIME(II) = TABLE(I)%X(1)%VALUES(II)*XGSCALE
450 YI(II) = TABLE(I)%Y%VALUES(II)*YGSCALE
451 ENDDO
452 ! Normalizing the curves to improve the optimization
453 MXVALT = MAXVAL(TIME)
454 MXVALY = MAXVAL(YI)
455 ! Considering infinite value shape
456 IF (ISHAPE == 1) THEN
457 GINFINI = YI(SIZE_TIME)/MXVALY
458 ELSE
459 GINFINI = ZERO
460 ENDIF
461 TIME = TIME/MXVALT
462 YI = YI/MXVALY
463 FOUND = .TRUE.
464 EXIT
465 ENDIF
466 ENDDO
467.NOT. IF (FOUND) THEN
468 CALL ANCMSG(MSGID=1928,MSGTYPE=MSGERROR,
469 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,I2=FCT_ID)
470 ENDIF
471C
472 ! Order of prony series to high for the number of data points
473 IF (2*NPRONY >= SIZE_TIME) THEN
474 CALL ANCMSG(MSGID=1921,MSGTYPE=MSGERROR,
475 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,I2=NPRONY,
476 . I3=SIZE_TIME,I4=FLOOR(SIZE_TIME/TWO))
477 ENDIF
478C
479 ! Allocation and initialization of the vectors
480 ALLOCATE(XI(2*NPRONY),XI_TR(2*NPRONY),JAC(SIZE_TIME,2*NPRONY),FXI(SIZE_TIME),
481 . DELTAY(SIZE_TIME),JACTJAC(2*NPRONY,2*NPRONY),ID(2*NPRONY,2*NPRONY),
482 . JACTDY(2*NPRONY),IPIV(2*NPRONY),JACT(2*NPRONY,SIZE_TIME),FXI_TR(SIZE_TIME),
483 . VECT(2*NPRONY),MU(2*NPRONY))
484C
485 ! Initialization of parameters value
486 DO I = 1,NPRONY
487 XI(2*I-1) = ONE/(I+ISHAPE)
488 XI(2*I) = I-1+ISHAPE
489 ENDDO
490C
491 ! Identity matrix
492 ID = ZERO
493 DO I = 1,2*NPRONY
494 ID(I,I) = ONE
495 ENDDO
496C
497 ! Initialization of the residues
498 FXI(:) = GINFINI
499 DO I = 1,SIZE_TIME
500 DO II = 1,NPRONY
501 FXI(I) = FXI(I) + XI(2*II-1)*EXP(-XI(2*II)*TIME(I))
502 JAC(I,2*II-1) = -EXP(-XI(2*II)*TIME(I))
503 JAC(I,2*II) = TIME(I)*XI(2*II-1)*EXP(-XI(2*II)*TIME(I))
504 ENDDO
505 ENDDO
506 DELTAY = YI - FXI
507C
508 ! Computation of the squared sum of residues
509 RI2 = DOT_PRODUCT(DELTAY,DELTAY)
510C
511 ! Initialization of the damping parameter and the convergence criterion
512 LAMBDA = MAXVAL(JAC)
513 NU = TWO
514 CONV = .FALSE.
515 NITER = 0
516C
517 !-------------------------------------------------------------------------------------
518 ! Levenberg Marquardt Algorithm + Line-search
519 !-------------------------------------------------------------------------------------
520#ifndef WITHOUT_LINALG
521.NOT..AND. DO WHILE ((CONV)(NITER<MXITER))
522C
523 ! a) Compute the Jacobian matrix products
524 ! -> Transpose of Jacobian matrix
525 JACT = TRANSPOSE(JAC)
526 ! -> JtJ product
527 JACTJAC = MATMUL(JACT,JAC)
528 ! -> Adding the damping parameter
529 JACTJAC = JACTJAC + LAMBDA*ID
530 ! -> Compute the gradient of residues
531 JACTDY = -MATMUL(JACT,DELTAY)
532c
533 ! b) Compute the trial values of XI
534 CALL DGESV(2*NPRONY, 1, JACTJAC, 2*NPRONY, IPIV, JACTDY , 2*NPRONY, INFO)
535 XI_TR = XI + JACTDY
536C
537 ! Checking convergence
538 IF (SQRT(DOT_PRODUCT(JACTDY,JACTDY)) <= EPS2*(SQRT(DOT_PRODUCT(XI,XI))+EPS2)) THEN
539 CONV = .TRUE.
540 ELSE
541 ! c) Compute the trial values of the function
542 FXI_TR(:) = GINFINI
543 DO I = 1,SIZE_TIME
544 DO II = 1,NPRONY
545 FXI_TR(I) = FXI_TR(I) + XI_TR(2*II-1)*EXP(-XI_TR(2*II)*TIME(I))
546 ENDDO
547 ENDDO
548 ! d) Computation of the squared sum of residues
549 RI2_OLD = RI2
550 RI2 = DOT_PRODUCT(YI-FXI_TR,YI-FXI_TR)
551C
552 ! Compute gain ration
553 VECT = MATMUL(JACT,DELTAY)
554 VECT = LAMBDA*JACTDY - VECT
555 DENOM = DOT_PRODUCT(JACTDY,VECT)
556 RHO = (RI2_OLD - RI2) / MAX(DENOM,EM20)
557C
558 ! e) Update parameters
559 ! -> Step accepted
560.AND. IF ((RHO > ZERO)(INFO == 0)) THEN
561 ! Update damping parameter
562 LAMBDA = LAMBDA*MAX(THIRD,ONE-(TWO*RHO-ONE)**3)
563 NU = TWO
564 ! Update XI vector with line search to be sure that XI >= ZERO
565 ! -> Compute the MU vector
566 DO I = 1,2*NPRONY
567 IF (JACTDY(I)>=ZERO) THEN
568 MU(I) = ONE
569 ELSE
570 MU(I) = MIN(ONE,((ONE-ALPHA)/(-JACTDY(I)))*XI(I))
571 ENDIF
572 XI(I) = XI(I) + MU(I)*JACTDY(I)
573 ENDDO
574 ! Update FXI vector
575 FXI(:) = GINFINI
576 DO I = 1,SIZE_TIME
577 DO II = 1,NPRONY
578 FXI(I) = FXI(I) + XI(2*II-1)*EXP(-XI(2*II)*TIME(I))
579 ENDDO
580 ENDDO
581 ! Update residue vector
582 DELTAY = YI - FXI
583 RI2 = DOT_PRODUCT(DELTAY,DELTAY)
584 ! Update the Jacobian matrix
585 JAC = ZERO
586 DO I = 1,SIZE_TIME
587 DO II = 1,NPRONY
588 JAC(I,2*II-1) = -EXP(-XI(2*II)*TIME(I))
589 JAC(I,2*II) = TIME(I)*XI(2*II-1)*EXP(-XI(2*II)*TIME(I))
590 ENDDO
591 ENDDO
592 JACT = TRANSPOSE(JAC)
593 ! -> Step not accepted
594 ELSE
595 ! Update damping parameter
596 LAMBDA = LAMBDA*NU
597 NU = NU*TWO
598 ENDIF
599C
600 ! f) Computing the convergence criterion
601 JACTDY = TWO*MATMUL(JACT,DELTAY)
602 JACTDY_N2 = SQRT(DOT_PRODUCT(JACTDY,JACTDY))
603 IF (JACTDY_N2 < EPS1) CONV = .TRUE.
604c
605 ! g) Increasing the iteration number
606 NITER = NITER + 1
607c
608 ENDIF
609C
610 ENDDO
611#else
612 CONV = .FALSE.
613 WRITE(6,*) "Error: Blas/Lapack required"
614#endif
615
616C
617 ! Algorithm has not converged
618.OR..NOT. IF (NITER == MXITER CONV) THEN
619 CALL ANCMSG(MSGID=1926,MSGTYPE=MSGERROR,
620 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID)
621 ! Solutions values must be checked
622.OR..OR. ELSEIF (ABS(JACTDY_N2)>EPS1 MAXVAL(XI)>EP20 MINVAL(XI)<EM20) THEN
623 CALL ANCMSG(MSGID=1927,MSGTYPE=MSGWARNING,
624 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID)
625 ENDIF
626C
627 ! Storage of the results in Prony parameters (accounting for normalizing values)
628 G = ZERO
629 BETA = ZERO
630 IF (ISHAPE == 1) THEN
631 GINFINI = GINFINI*MXVALY
632 G(1) = GINFINI
633 BETA(1) = ZERO
634 ENDIF
635 DO I = 1,NPRONY
636 G(I+ISHAPE) = XI(2*I-1)*MXVALY
637 BETA(I+ISHAPE) = XI(2*I)/MXVALT
638 ENDDO
639 COST = RI2
640 DERIV = JACTDY_N2
641C
642 ! Deallocation of tables
643 DEALLOCATE(XI,XI_TR,JAC,FXI,DELTAY,JACTJAC,ID,JACTDY,IPIV,JACT,FXI_TR,VECT,YI,TIME,MU)
644 END
645
646!||====================================================================
647!|| lm_least_square_prony_2 ../starter/source/materials/visc/hm_read_visc_prony.F
648!||--- called by ------------------------------------------------------
649!|| hm_read_visc_prony ../starter/source/materials/visc/hm_read_visc_prony.F
650!||--- calls -----------------------------------------------------
651!|| ancmsg ../starter/source/output/message/message.F
652!|| qsort ../starter/source/materials/visc/hm_read_visc_prony.F
653!||--- uses -----------------------------------------------------
654!|| message_mod ../starter/share/message_module/message_mod.F
655!|| table_mod ../starter/share/modules1/table_mod.F
656!||====================================================================
657 SUBROUTINE LM_LEAST_SQUARE_PRONY_2(MAT_ID ,NPRONY ,TABLE ,FCT_IDS ,XGS_SCALE,YGS_SCALE,
658 . FCT_IDL ,XGL_SCALE,YGL_SCALE,G ,BETA ,COST ,
659 . DERIV ,ISHAPE ,GINFINI )
660C-----------------------------------------------
661C M o d u l e s
662C-----------------------------------------------
663 USE TABLE_MOD
664 USE MESSAGE_MOD
665C-----------------------------------------------
666C I m p l i c i t T y p e s
667C-----------------------------------------------
668#include "implicit_f.inc"
669#include "com04_c.inc"
670C-----------------------------------------------
671C Dummy arguments
672C-----------------------------------------------
673 INTEGER, INTENT(IN) :: NPRONY,FCT_IDS,FCT_IDL,MAT_ID,ISHAPE
674 my_real
675 . XGS_SCALE,YGS_SCALE,XGL_SCALE,YGL_SCALE,COST,DERIV,GINFINI
676 my_real, DIMENSION(100) :: G,BETA
677 TYPE(TTABLE) TABLE(NTABLE)
678C-----------------------------------------------
679C Local arguments
680C-----------------------------------------------
681 INTEGER I,II,INFO,SIZE_FREQ,NITER,ITABS,ITABL,FLAG,
682 . SIZE_FREQS,SIZE_FREQL,K
683 INTEGER, DIMENSION(:), ALLOCATABLE :: IPIV,IDX
684 DOUBLE PRECISION ::
685 . RI2,RI2_OLD,LAMBDA,JACTDY_N2,
686 . RHO,NU,DENOM,XI_N2,MXFREQ,MXNORM
687 DOUBLE PRECISION, PARAMETER :: EPS1 = 1.0D-8
688 DOUBLE PRECISION, PARAMETER :: EPS2 = 1.0D-10
689 DOUBLE PRECISION, PARAMETER :: ALPHA = 0.5D0
690 DOUBLE PRECISION, PARAMETER :: MXITER = 10000
691 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: XI,FXIS,FXIL,DELTAYS,DELTAYL,
692 . JACTDY,YIS,YIL,FREQ,XI_TR,FXIS_TR,
693 . FXIL_TR,VECT,MU,FREQS,FREQL,YIS_TEMP,
694 . YIL_TEMP
695 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: JACS,JACL,JACTJAC,JACST,JACLT,ID,WEIGHT
696 LOGICAL :: CONV
697C
698 ! Finding the functions for Storage and Loss moduli
699 ITABS = 0
700 ITABL = 0
701 DO I = 1,NTABLE
702 IF (TABLE(I)%NOTABLE == FCT_IDS) ITABS = I
703 IF (TABLE(I)%NOTABLE == FCT_IDL) ITABL = I
704.AND. IF (ITABS /= 0 ITABL /= 0) EXIT
705 ENDDO
706 IF (ITABS == 0) THEN
707 CALL ANCMSG(MSGID=1924,MSGTYPE=MSGERROR,
708 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,I2=FCT_IDS)
709 ENDIF
710 IF (ITABL == 0) THEN
711 CALL ANCMSG(MSGID=1925,MSGTYPE=MSGERROR,
712 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,I2=FCT_IDL)
713 ENDIF
714C
715 ! --------------------------------------------------------------------------
716 ! Recover the experimental data and unifying the abscissa vector if needed
717 ! --------------------------------------------------------------------------
718 ! -> Temporary data storage
719 ! --------------------------------------------------------------------
720 SIZE_FREQS = SIZE(TABLE(ITABS)%X(1)%VALUES)
721 SIZE_FREQL = SIZE(TABLE(ITABL)%X(1)%VALUES)
722 ALLOCATE(FREQS(SIZE_FREQS),YIS_TEMP(SIZE_FREQS),
723 . FREQL(SIZE_FREQL),YIL_TEMP(SIZE_FREQL))
724 FREQS(1:SIZE_FREQS) = TABLE(ITABS)%X(1)%VALUES(1:SIZE_FREQS)*XGS_SCALE
725 YIS_TEMP(1:SIZE_FREQS) = TABLE(ITABS)%Y%VALUES(1:SIZE_FREQS)*YGS_SCALE
726 FREQL(1:SIZE_FREQL) = TABLE(ITABL)%X(1)%VALUES(1:SIZE_FREQL)*XGL_SCALE
727 YIL_TEMP(1:SIZE_FREQL) = TABLE(ITABL)%Y%VALUES(1:SIZE_FREQL)*YGL_SCALE
728C
729 ! -> Number of abscissa is the same
730 ! --------------------------------------------------------------------
731 FLAG = 0
732 SIZE_FREQ = 0
733 IF (SIZE_FREQS == SIZE_FREQL) THEN
734 ! If there is no difference at all between abscissa
735 IF (SQRT(DOT_PRODUCT(FREQS-FREQL,FREQS-FREQL))<EM08) THEN
736 FLAG = 0
737 SIZE_FREQ = SIZE_FREQS
738 ALLOCATE(FREQ(SIZE_FREQ),YIS(SIZE_FREQ),YIL(SIZE_FREQ))
739 MXFREQ = MAXVAL(FREQS)
740 DO II = 1, SIZE_FREQ
741 FREQ(II) = FREQS(II)*(HUNDRED/MXFREQ)
742 YIS(II) = YIS_TEMP(II)
743 YIL(II) = YIL_TEMP(II)
744 ENDDO
745 ! If they have the same number but different value
746 ELSE
747 FLAG = 1
748 ENDIF
749 ENDIF
750 ! -> Abscissa are not the same
751 ! --------------------------------------------------------------------
752.OR. IF (SIZE_FREQ /= SIZE_FREQL FLAG == 1) THEN
753 ! -> Counting the total number of abscissa without duplicate identical values
754 !----------------------------------------------------------------------------
755 SIZE_FREQ = SIZE_FREQS + SIZE_FREQL
756 DO I = 1,SIZE_FREQL
757 FLAG = 0
758 DO II = 1,SIZE_FREQS
759 IF (ABS(FREQS(II)-FREQL(I))<EM08) THEN
760 FLAG = 1
761 ENDIF
762 ENDDO
763 IF (FLAG == 1) SIZE_FREQ = SIZE_FREQ - 1
764 ENDDO
765 ! -> Allocating and filling tables
766 !----------------------------------------------------------------------------
767 ALLOCATE(IDX(SIZE_FREQ),FREQ(SIZE_FREQ),YIS(SIZE_FREQ),YIL(SIZE_FREQ))
768 FREQ(1:SIZE_FREQS) = FREQS(1:SIZE_FREQS)
769 K = 0
770 DO I = 1,SIZE_FREQL
771 FLAG = 0
772 DO II = 1,SIZE_FREQS
773 IF (ABS(FREQS(II)-FREQL(I))<EM08) THEN
774 FLAG = 1
775 K = K + 1
776 ENDIF
777 ENDDO
778 IF (FLAG == 0) FREQ(SIZE_FREQS+I-K) = FREQL(I)
779 ENDDO
780 ! -> Sorting abscissa
781 !--------------------
782 DO II = 1,SIZE_FREQ
783 IDX(II) = II
784 ENDDO
785 CALL QSORT(FREQ,IDX,1,SIZE_FREQ)
786 MXFREQ = MAXVAL(FREQ)
787 FREQ = FREQ*(HUNDRED/MXFREQ)
788 FREQS = FREQS*(HUNDRED/MXFREQ)
789 FREQL = FREQL*(HUNDRED/MXFREQ)
790 ! -> Data interpolation with new abscissa
791 !----------------------------------------------------------------------------
792 II = 1
793 DO I = 1,SIZE_FREQ
794 IF (FREQ(I)<FREQS(1)) THEN
795 YIS(I) = YIS_TEMP(1) + ((YIS_TEMP(2)-YIS_TEMP(1))/(FREQS(2)-FREQS(1)))*(FREQ(I)-FREQS(1))
796 ELSEIF (FREQ(I)>FREQS(SIZE_FREQS)) THEN
797 YIS(I) = YIS_TEMP(SIZE_FREQS-1) + ((YIS_TEMP(SIZE_FREQS)-YIS_TEMP(SIZE_FREQS-1))/
798 . (FREQS(SIZE_FREQS)-FREQS(SIZE_FREQS-1)))*(FREQ(I)-FREQS(SIZE_FREQS-1))
799 ELSE
800 IF (FREQ(I)>FREQS(II+1)) THEN
801 II = II + 1
802 II = MIN(II,SIZE_FREQS-1)
803 ENDIF
804 YIS(I) = YIS_TEMP(II) + ((YIS_TEMP(II+1)-YIS_TEMP(II))/(FREQS(II+1)-FREQS(II)))*(FREQ(I)-FREQS(II))
805 ENDIF
806 ENDDO
807 II = 1
808 DO I = 1,SIZE_FREQ
809 IF (FREQ(I)<FREQL(1)) THEN
810 YIL(I) = YIL_TEMP(1) + ((YIL_TEMP(2)-YIL_TEMP(1))/(FREQL(2)-FREQL(1)))*(FREQ(I)-FREQL(1))
811 ELSEIF (FREQ(I)>FREQL(SIZE_FREQL)) THEN
812 YIL(I) = YIL_TEMP(SIZE_FREQL-1) + ((YIL_TEMP(SIZE_FREQL)-YIL_TEMP(SIZE_FREQL-1))/
813 . (FREQL(SIZE_FREQL)-FREQL(SIZE_FREQL-1)))*(FREQ(I)-FREQL(SIZE_FREQL-1))
814 ELSE
815 IF (FREQ(I)>FREQL(II+1)) THEN
816 II = II + 1
817 II = MIN(II,SIZE_FREQL-1)
818 ENDIF
819 YIL(I) = YIL_TEMP(II) + ((YIL_TEMP(II+1)-YIL_TEMP(II))/(FREQL(II+1)-FREQL(II)))*(FREQ(I)-FREQL(II))
820 ENDIF
821 ENDDO
822 ENDIF
823C
824 ! Order of prony series to high for the number of data points
825 IF (2*NPRONY >= SIZE_FREQ) THEN
826 CALL ANCMSG(MSGID=1921,MSGTYPE=MSGERROR,
827 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,I2=NPRONY,
828 . I3=SIZE_FREQ,I4=INT(SIZE_FREQ/TWO))
829 ENDIF
830C
831 ! Allocation and initialization of the vectors
832 ALLOCATE(XI(2*NPRONY),XI_TR(2*NPRONY),JACS(SIZE_FREQ,2*NPRONY),JACL(SIZE_FREQ,2*NPRONY),
833 . FXIS(SIZE_FREQ),FXIL(SIZE_FREQ),DELTAYS(SIZE_FREQ),DELTAYL(SIZE_FREQ),
834 . JACTJAC(2*NPRONY,2*NPRONY),ID(2*NPRONY,2*NPRONY),JACTDY(2*NPRONY),IPIV(2*NPRONY),
835 . JACST(2*NPRONY,SIZE_FREQ),JACLT(2*NPRONY,SIZE_FREQ),FXIS_TR(SIZE_FREQ),
836 . FXIL_TR(SIZE_FREQ),VECT(2*NPRONY),MU(2*NPRONY),WEIGHT(SIZE_FREQ,SIZE_FREQ))
837C
838 ! Initialization of parameters value
839 DO I = 1,NPRONY
840 XI(2*I-1) = YIS(1)/(I+ISHAPE)
841 XI(2*I) = I-1+ISHAPE
842 ENDDO
843C
844 ! Identity matrix
845 ID = ZERO
846 DO I = 1,2*NPRONY
847 ID(I,I) = ONE
848 ENDDO
849C
850 ! Considering infinite value shape
851 IF (ISHAPE == 1) THEN
852 GINFINI = YIS(1)
853 ELSE
854 GINFINI = ZERO
855 ENDIF
856C
857 ! Weight matrix
858 WEIGHT = ZERO
859 MXNORM = MAX(MAXVAL(YIS),MAXVAL(YIL))
860 DO I = 1,SIZE_FREQ
861 WEIGHT(I,I) = ONE/(MXNORM**2)
862 ENDDO
863C
864 ! Initialization of the Prony series and the residues
865 FXIS(:) = GINFINI
866 FXIL(:) = ZERO
867 DO I = 1,SIZE_FREQ
868 DO II = 1,NPRONY
869 ! Storage modulus series
870 FXIS(I) = FXIS(I) + ((XI(2*II-1)*(FREQ(I)**2))/(XI(2*II)**2 + FREQ(I)**2))
871 JACS(I,2*II-1) = -(FREQ(I)**2)/(XI(2*II)**2 + FREQ(I)**2)
872 JACS(I,2*II) = TWO*XI(2*II-1)*XI(2*II)*(FREQ(I)**2)/((XI(2*II)**2 + FREQ(I)**2)**2)
873 ! Loss modulus series
874 FXIL(I) = FXIL(I) + ((XI(2*II-1)*XI(2*II)*FREQ(I))/(XI(2*II)**2 + FREQ(I)**2))
875 JACL(I,2*II-1) = -(FREQ(I)*XI(2*II))/(XI(2*II)**2 + FREQ(I)**2)
876 JACL(I,2*II) = XI(2*II-1)*FREQ(I)*((XI(2*II)**2 - FREQ(I)**2)/(XI(2*II)**2 + FREQ(I)**2)**2)
877 ENDDO
878 ENDDO
879 DELTAYS = YIS - FXIS
880 DELTAYL = YIL - FXIL
881C
882 ! Computation of the squared sum of residues
883 RI2 = ((ONE/MXNORM)**2)*(DOT_PRODUCT(DELTAYS,DELTAYS) + DOT_PRODUCT(DELTAYL,DELTAYL))
884C
885 ! Initialization of the damping parameter and the convergence criterion
886 LAMBDA = MAXVAL(JACS) + MAXVAL(JACL)
887 NU = TWO
888 CONV = .FALSE.
889 NITER = 0
890 JACTDY_N2 = ZERO
891C
892 !-------------------------------------------------------------------------------------
893 ! Levenberg Marquardt Algorithm + Line-search
894 !-------------------------------------------------------------------------------------
895#ifndef WITHOUT_LINALG
896.NOT..AND. DO WHILE ((CONV)(NITER<MXITER))
897C
898 ! a) Compute the Jacobian matrix products
899 ! -> Transpose of Jacobian matrix
900 JACST = TRANSPOSE(JACS)
901 JACLT = TRANSPOSE(JACL)
902 ! -> JtJ product
903 JACTJAC = MATMUL(JACST,MATMUL(WEIGHT,JACS))
904 JACTJAC = JACTJAC + MATMUL(JACLT,MATMUL(WEIGHT,JACL))
905 ! -> Adding the damping parameter
906 JACTJAC = JACTJAC + LAMBDA*ID
907 ! -> Compute the gradient of residues
908 JACTDY = -MATMUL(JACST,MATMUL(WEIGHT,DELTAYS))
909 JACTDY = JACTDY - MATMUL(JACLT,MATMUL(WEIGHT,DELTAYL))
910c
911 ! b) Compute the trial values of XI
912 CALL DGESV(2*NPRONY, 1, JACTJAC, 2*NPRONY, IPIV, JACTDY , 2*NPRONY, INFO)
913 XI_TR = XI + JACTDY
914C
915 ! Checking convergence
916 IF (SQRT(DOT_PRODUCT(JACTDY,JACTDY)) <= EPS2*(SQRT(DOT_PRODUCT(XI,XI))+EPS2)) THEN
917 CONV = .TRUE.
918 ELSE
919 ! c) Compute the trial values of the function
920 FXIS_TR(:) = GINFINI
921 FXIL_TR(:) = ZERO
922 DO I = 1,SIZE_FREQ
923 DO II = 1,NPRONY
924 FXIS_TR(I) = FXIS_TR(I) + ((XI_TR(2*II-1)*(FREQ(I)**2))/(XI_TR(2*II)**2 + FREQ(I)**2))
925 FXIL_TR(I) = FXIL_TR(I) + ((XI_TR(2*II-1)*XI_TR(2*II)*FREQ(I))/(XI_TR(2*II)**2 + FREQ(I)**2))
926 ENDDO
927 ENDDO
928 ! d) Computation of the squared sum of residues
929 RI2_OLD = RI2
930 RI2 = ((ONE/MXNORM)**2)*(DOT_PRODUCT((YIS-FXIS_TR),(YIS-FXIS_TR)) +
931 . DOT_PRODUCT((YIL-FXIL_TR),(YIL-FXIL_TR)))
932C
933 ! Compute gain ratio
934 VECT = MATMUL(JACST,MATMUL(WEIGHT,DELTAYS)) + MATMUL(JACLT,MATMUL(WEIGHT,DELTAYL))
935 VECT = LAMBDA*JACTDY - VECT
936 DENOM = DOT_PRODUCT(JACTDY,VECT)
937 RHO = (RI2_OLD - RI2) / MAX(DENOM,EM20)
938C
939 ! e) Update parameters
940 ! -> Step accepted
941.AND. IF ((RHO > ZERO)(INFO == 0)) THEN
942 ! Update damping parameter
943 LAMBDA = LAMBDA*MAX(THIRD,ONE-(TWO*RHO-ONE)**3)
944 NU = TWO
945 ! Update XI vector
946 ! -> Compute the MU vector for line search
947 DO I = 1,2*NPRONY
948 IF (JACTDY(I)>=ZERO) THEN
949 MU(I) = ONE
950 ELSE
951 MU(I) = MIN(ONE,((ONE-ALPHA)/(-JACTDY(I)))*XI(I))
952 ENDIF
953 XI(I) = XI(I) + MU(I)*JACTDY(I)
954 ENDDO
955 ! Update FXI vector
956 FXIS(:) = GINFINI
957 FXIL(:) = ZERO
958 DO I = 1,SIZE_FREQ
959 DO II = 1,NPRONY
960 FXIS(I) = FXIS(I) + ((XI(2*II-1)*(FREQ(I)**2))/(XI(2*II)**2 + FREQ(I)**2))
961 FXIL(I) = FXIL(I) + ((XI(2*II-1)*XI(2*II)*FREQ(I))/(XI(2*II)**2 + FREQ(I)**2))
962 ENDDO
963 ENDDO
964 ! Update residue vector
965 DELTAYS = YIS - FXIS
966 DELTAYL = YIL - FXIL
967 RI2 = ((ONE/MXNORM)**2)*(DOT_PRODUCT(DELTAYS,DELTAYS) + DOT_PRODUCT(DELTAYL,DELTAYL))
968 ! Update the Jacobian matrix
969 JACS = ZERO
970 DO I = 1,SIZE_FREQ
971 DO II = 1,NPRONY
972 JACS(I,2*II-1) = -(FREQ(I)**2)/(XI(2*II)**2 + FREQ(I)**2)
973 JACS(I,2*II) = TWO*XI(2*II-1)*XI(2*II)*(FREQ(I)**2)/((XI(2*II)**2 + FREQ(I)**2)**2)
974 JACL(I,2*II-1) = -(FREQ(I)*XI(2*II))/(XI(2*II)**2 + FREQ(I)**2)
975 JACL(I,2*II) = XI(2*II-1)*FREQ(I)*((XI(2*II)**2 - FREQ(I)**2)/(XI(2*II)**2 + FREQ(I)**2)**2)
976 ENDDO
977 ENDDO
978 JACST = TRANSPOSE(JACS)
979 JACLT = TRANSPOSE(JACL)
980 ! -> Step not accepted
981 ELSE
982 ! Update damping parameter
983 LAMBDA = LAMBDA*NU
984 NU = NU*TWO
985 ENDIF
986C
987 ! f) Computing the convergence criterion
988 JACTDY = TWO*MATMUL(JACST,MATMUL(WEIGHT,DELTAYS)) + TWO*MATMUL(JACLT,MATMUL(WEIGHT,DELTAYL))
989 JACTDY_N2 = SQRT(DOT_PRODUCT(JACTDY,JACTDY))
990 IF (JACTDY_N2 < EPS1) CONV = .TRUE.
991c
992 ! g) Increasing the iteration number
993 NITER = NITER + 1
994c
995 ENDIF
996C
997 ENDDO
998#else
999 CONV = .FALSE.
1000 WRITE(6,*) "Error: Blas/Lapack required"
1001#endif
1002
1003C
1004 ! Algorithm has not converged
1005.OR..NOT. IF (NITER == MXITER CONV) THEN
1006 CALL ANCMSG(MSGID=1926,MSGTYPE=MSGERROR,
1007 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID)
1008 ! Solutions values must be checked
1009.OR..OR. ELSEIF (ABS(JACTDY_N2)>EPS1 MAXVAL(XI)>EP10 MINVAL(XI)<EM10) THEN
1010 CALL ANCMSG(MSGID=1927,MSGTYPE=MSGWARNING,
1011 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID)
1012 ENDIF
1013C
1014 ! Storage of the results in Prony parameters (accounting for normalizing values)
1015 G = ZERO
1016 BETA = ZERO
1017 IF (ISHAPE == 1) THEN
1018 G(1) = GINFINI
1019 BETA(1) = ZERO
1020 ENDIF
1021 DO I = 1,NPRONY
1022 G(I+ISHAPE) = XI(2*I-1)
1023 BETA(I+ISHAPE) = XI(2*I)*(MXFREQ/HUNDRED)
1024 ENDDO
1025 COST = RI2
1026 DERIV = JACTDY_N2
1027C
1028 ! Deallocation of tables
1029 DEALLOCATE(XI,XI_TR,JACS,JACL,FXIS,FXIL,DELTAYS,DELTAYL,JACTJAC,ID,JACTDY,IPIV,
1030 . JACST,JACLT,FXIS_TR,FXIL_TR,VECT,MU,FREQS,YIS_TEMP,FREQL,YIL_TEMP)
1031 IF (ALLOCATED(IDX)) DEALLOCATE(IDX)
1032 END
1033
1034
1035C---------------------------------------------------------------------------------------
1036!||====================================================================
1037!|| qsort ../starter/source/materials/visc/hm_read_visc_prony.F
1038!||--- called by ------------------------------------------------------
1039!|| lm_least_square_prony_2 ../starter/source/materials/visc/hm_read_visc_prony.F
1040!||--- calls -----------------------------------------------------
1041!||====================================================================
1042 RECURSIVE SUBROUTINE QSORT(A, IDX, FIRST, LAST)
1043C-----------------------------------------------
1044C I m p l i c i t T y p e s
1045C-----------------------------------------------
1046#include "implicit_f.inc"
1047C-----------------------------------------------
1048C D u m m y A r g u m e n t s
1049C-----------------------------------------------
1050 DOUBLE PRECISION, INTENT(INOUT) :: A(*)
1051 INTEGER, INTENT(IN) :: FIRST, LAST
1052 INTEGER, INTENT(INOUT) :: IDX(*)
1053C-----------------------------------------------
1054C L o c a l V a r i a b l e s
1055C-----------------------------------------------
1056 DOUBLE PRECISION ::
1057 . X, T
1058 INTEGER :: I, J, I1, I2
1059C-----------------------------------------------
1060C P r e - C o n d i t i o n
1061C-----------------------------------------------
1062 IF(FIRST>LAST)RETURN
1063C-----------------------------------------------
1064C S o u r c e L i n e s
1065C-----------------------------------------------
1066 X = A( (FIRST + LAST) / 2 )
1067 I = FIRST
1068 J = LAST
1069 DO
1070 DO WHILE (A(I) < X)
1071 I = I + 1
1072 ENDDO
1073 DO WHILE(X < A(J))
1074 J = J - 1
1075 ENDDO
1076 IF (I >= J) EXIT
1077 T = A(I)
1078 A(I) = A(J)
1079 A(J) = T
1080 I1 = IDX(i)
1081 IDX(I) = IDX(J)
1082 IDX(J) = I1
1083 I = I + 1
1084 J = J - 1
1085 ENDDO
1086 IF (FIRST < I - 1) CALL QSORT(A, IDX, FIRST, I - 1)
1087 IF (J + 1 < LAST) CALL QSORT(A, IDX, J + 1, LAST)
1088 END SUBROUTINE
my_real
#define my_real
Definition cppsort.cpp:32
hm_get_float_array_index
subroutine hm_get_float_array_index(name, rval, index, is_available, lsubmodel, unitab)
Definition hm_get_float_array_index.F:105
hm_get_floatv
subroutine hm_get_floatv(name, rval, is_available, lsubmodel, unitab)
Definition hm_get_floatv.F:415
hm_get_floatv_dim
subroutine hm_get_floatv_dim(name, dim_fac, is_available, lsubmodel, unitab)
Definition hm_get_floatv_dim.F:159
hm_get_intv
subroutine hm_get_intv(name, ival, is_available, lsubmodel)
Definition hm_get_intv.F:523
hm_option_is_encrypted
subroutine hm_option_is_encrypted(is_encrypted)
Definition hm_option_is_encrypted.F:284
hm_read_visc_prony
subroutine hm_read_visc_prony(visc, ivisc, ntable, table, mat_id, unitab, lsubmodel)
Definition hm_read_visc_prony.F:44
lm_least_square_prony
subroutine lm_least_square_prony(mat_id, nprony, table, fct_id, xgscale, ygscale, g, beta, cost, deriv, ishape, ginfini)
Definition hm_read_visc_prony.F:404
lm_least_square_prony_2
subroutine lm_least_square_prony_2(mat_id, nprony, table, fct_ids, xgs_scale, ygs_scale, fct_idl, xgl_scale, ygl_scale, g, beta, cost, deriv, ishape, ginfini)
Definition hm_read_visc_prony.F:660
id
initmumps id
Definition multiplerhs_example.m:4
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
hm_option_read_mod
Definition hm_option_read_mod.F:408
message_mod
Definition message_mod.F:1249
submodel_mod
Definition submodel_mod.F:712
table_mod
Definition table_mod.F:112
unitab_mod
Definition unitab_mod.F:513
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
starter
program starter
Definition starter.F:39
tabulated
subroutine tabulated(iflag, nel, pm, off, eint, mu, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, npf, tf)
Definition tabulated.F:32