OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
law69_nlsqf_auto.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!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
25!||--- called by ------------------------------------------------------
26!|| law69_upd ../starter/source/materials/mat/mat069/law69_upd.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| law69_check ../starter/source/materials/tools/nlsqf.F
30!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
31!|| law69_guess_bounds ../starter/source/materials/tools/nlsqf.f
32!|| mrqmin ../starter/source/materials/tools/nlsqf.F
33!|| ogden0 ../starter/source/materials/tools/nlsqf.F
34!||--- uses -----------------------------------------------------
35!|| message_mod ../starter/share/message_module/message_mod.F
36!||====================================================================
37 SUBROUTINE law69_nlsqf_auto(LAWID,STRETCH,Y,MA,NMUAL,NPT,MUAL,
38 $ ICHECK, NSTART, ERRTOL,ID,TITR , IS_ENCRYPTED)
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
42 USE message_mod
43 USE names_and_titles_mod , ONLY : nchartitle
44 !-----------------------------------------------
45 ! I m p l i c i t T y p e s
46 !-----------------------------------------------
47#include "implicit_f.inc"
48 !-----------------------------------------------
49 ! C o m m o n B l o c k s
50 !-----------------------------------------------
51#include "units_c.inc"
52 !-----------------------------------------------
53 ! D u m m y A r g u m e n t s
54 !-----------------------------------------------
55 INTEGER , INTENT(INOUT) :: LAWID, NMUAL, MA,ICHECK
56 INTEGER , INTENT(IN) :: NPT, NSTART,ID,IS_ENCRYPTED
57 my_real , INTENT(IN) :: errtol
58 my_real , DIMENSION(NPT) , INTENT(IN) :: y,stretch
59 my_real , DIMENSION(10), INTENT(INOUT) :: mual
60 CHARACTER(LEN=NCHARTITLE),INTENT(IN) :: TITR
61 !-----------------------------------------------
62 ! L o c a l V a r i a b l e s
63 !-----------------------------------------------
64 INTEGER MAXA
65 parameter(maxa=20)
66
67 INTEGER I,IDUM,ITER,ICOUNT,J,K,IRET ,ISTART, NPSAVED,
68 . ivalid,icomp, mmax,miniter_lm, maxiter_lm,
69 . cnt_hit_eps_lm,lmt_hit_eps_lm, lmstop,ifuncs,
70 . nguess,nmuals,loop_istart ,icheck_guess,
71 . jcheck, ifit_success,iend_finding,iend_iter,
72 . mu_incr_guess,irnd1,iswitch
73 INTEGER NONZERO(MAXA)
74 my_real gamma,errnow,gasdev,errpre,errmin, errmax,
75 . errave, errave_min,yogd,xogd,eps_lm,max_abs_yi,
76 . small_fac_abs_yi, small_abs_yi,spready,gamma_stop,
77 . err
78 my_real a(maxa),covar(maxa,maxa),alpha(maxa,maxa),
79 . mcof_min(maxa), mcof_max(maxa),as(10),a0(maxa),
80 . atry(maxa)
81
82 my_real, DIMENSION(:), ALLOCATABLE :: sig
83!!=======================================================================
84 !
85 iswitch = 0
86 miniter_lm = 5
87 maxiter_lm = 100
88 eps_lm = em03
89 lmt_hit_eps_lm = 2
90 gamma_stop = ep15
91 !! if check the validity of initial guess (0=no, 1=yes)
92 !! we don't want to check, because invalid initial point
93 !! could converge to valid point
94 icheck_guess = 0
95 !! if enforce mu(i) < mu(i+1) in generating initial guess
96 !! we don't need this anymore since we are using random numbers
97 !! instead of loop through all
98 !! the combinations
99 mu_incr_guess =0
100
101 ALLOCATE (sig(1:npt))
102 icomp=0
103 IF(icheck < 0) icomp= 1
104 icheck = abs(icheck)
105 ! IF ABS(Y(I)) < SMALL_ABS_YI, use SMALL_ABS_YI to avoid
106 ! unnecessary numerical issues when avoid divided by small value
107 max_abs_yi = zero
108 DO i = 1, npt
109 max_abs_yi = max( max_abs_yi, abs(y(i)) )
110 ENDDO
111
112 small_fac_abs_yi = em3
113 small_abs_yi = max_abs_yi * small_fac_abs_yi
114 spready=em02
115 IF(icomp == 0) THEN
116 DO j=1,npt
117 sig(j)=spready*max(small_abs_yi, abs(y(j)) )
118 ENDDO
119 ELSE
120 DO j=1,npt
121 sig(j) = y(j)
122 IF(sig(j) == zero) sig(j) = em15
123 ENDDO
124 ENDIF
125C=======================================================================
126 npsaved = 0
127 as(1:10) = zero
128 ifit_success = 0
129 iend_iter = 0
130 iend_finding = 0
131 jcheck = 2 ! stating Criteria
132 ifuncs = 1
133 lmstop =0
134 errave_min = huge(errave_min)
135 DO WHILE(ifit_success == 0 ) ! f
136 IF(iend_iter <= 1) THEN
137 DO i=1,nmual
138 nonzero(i)=1
139 ENDDO
140 IF (lawid == 1) then ! OGDEN
141 DO i=nmual+1,10
142 nonzero(i)=0
143 a(i)=zero
144 ENDDO
145 ELSEIF (lawid == 2) then ! Mooney-Rivlin
146 DO i=nmual+1,ma
147 nonzero(i)=0
148 ENDDO
149 nonzero(2)=0
150 nonzero(4)=0
151 DO i=nmual+1,ma
152 a(i)=zero
153 ENDDO
154 ENDIF
155 ENDIF
156 !
157 istart = 0
158 !! calculate MCOF_MIN, MCOF_MAX
159 CALL law69_guess_bounds(lawid, jcheck, stretch, y, npt,
160 . nmual, mcof_min, mcof_max,icomp)
161
162 !! initialize random number generator in LAW69_GUESS1
163 irnd1 = 205
164 iend_iter = 0
165 atry(1:maxa) = zero ! Initialization of ATRY
166 DO WHILE ( istart < nstart ) !!! 111 continue
167 !! get one guess point A0(1:NMUAL)
168 CALL law69_guess1(
169 $ lawid, nmual, icheck_guess, mu_incr_guess,irnd1,
170 $ a0, nonzero, mcof_min, mcof_max, nguess)
171
172 IF (nguess /= 0) THEN
173C calculate averaged ERROR before LM
174 errave = zero
175 IF(icomp == 0) THEN
176 DO i=1,npt
177 CALL ogden0(stretch(i), a0, yogd, nmual)
178 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
179 errave = errave + err
180 ENDDO
181 ELSE
182 DO i=1,npt
183 CALL ogden0(stretch(i), a0, yogd, nmual)
184 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
185 errave = errave + err
186 ENDDO
187 ENDIF
188C
189 errave = errave / (one * npt)
190 istart = istart + 1
191 ! start loop of LM optimization
192 iter = 0
193 cnt_hit_eps_lm = 0
194
195 iter = iter + 1
196 gamma=-1
197 CALL mrqmin(stretch,y,sig,npt,a0,nonzero,
198 . nmual,covar,alpha,ma,errnow,
199 . ifuncs,gamma,iret,icomp,mmax,atry)
200 IF (iret /= 0) cycle
201 lmstop = 0
202 loop_istart = 0
203 DO WHILE (lmstop == 0 )
204 errpre=errnow
205 iter = iter + 1
206 CALL mrqmin(stretch,y,sig,npt,a0,nonzero,
207 . ma,covar,alpha, ma, errnow,
208 . ifuncs,gamma,iret,icomp,mmax,atry)
209 IF (iret /= 0)THEN
210 loop_istart = 1
211 EXIT
212 ENDIF
213 !! IF A0(J) is too small, restart from next initial guess
214 DO j = 1, nmual
215 IF ( abs( a0(j) ) < 1e-20 ) THEN
216 loop_istart = 1
217 EXIT
218 ENDIF
219 ENDDO
220 IF(loop_istart == 1) EXIT
221c check convergence of LM optimization
222 IF ( iter > miniter_lm ) THEN
223 IF (errnow < errpre) THEN
224 IF ( abs(errpre) > zero) THEN
225 IF ( abs( (errnow-errpre)/ errpre ) < eps_lm) THEN
226 cnt_hit_eps_lm = cnt_hit_eps_lm + 1
227 IF ( cnt_hit_eps_lm >= lmt_hit_eps_lm ) THEN
228 lmstop = 1
229 ENDIF
230 ENDIF
231 ENDIF
232 ELSEIF (iter >= maxiter_lm .OR. gamma >= gamma_stop) THEN
233 lmstop = 1
234 ENDIF
235 ENDIF
236 ENDDO ! LMSTOP
237 IF(loop_istart == 1) cycle
238!
239 errave = 0.0
240 IF(icomp == 0) THEN
241 DO i=1,npt
242 CALL ogden0(stretch(i), a0, yogd, nmual)
243 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
244 errave = errave + err
245 ENDDO
246 ELSE
247 DO i=1,npt
248 CALL ogden0(stretch(i), a0, yogd, nmual)
249 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
250 errave = errave + err
251 ENDDO
252 ENDIF
253 errave = errave / (one * npt)
254 CALL law69_check(lawid, a0, nmual, jcheck, 0, ivalid)
255 IF (ivalid > 0) THEN
256 IF (npsaved==0 .OR.
257 $ ( npsaved > 0 .AND. errave < errave_min) ) THEN
258 npsaved = npsaved + 1
259 errave_min = errave
260 nmuals = nmual
261 DO i = 1, nmual
262 a(i) = a0(i)
263 as(i) = a0(i)
264 ENDDO
265 ENDIF
266 ENDIF
267 ! Check the fit
268 IF (npsaved > 0 .AND. errave_min < errtol ) THEN
269 ifit_success = 1
270 EXIT
271 ENDIF
272 ENDIF ! NGUESS
273 iswitch = 0
274 IF(ifit_success == 0 .AND. istart < nstart ) THEN
275 cycle
276 ELSEIF(ifit_success == 0 .OR. (nguess == 0 .AND. ifit_success ==0)) THEN
277 iswitch = 1
278 IF (jcheck == 2) THEN
279 jcheck = 1
280 iend_iter = 2
281 ELSEIF(jcheck == 1) THEN
282 IF(lawid == 2) THEN
283 lawid = 1
284 jcheck = 2
285 iend_iter =1
286 ELSEIF(nmual < 10) THEN
287 nmual = min(10, nmual + 2)
288 jcheck = 2
289 iend_iter = 1
290 ELSEIF(icomp == 0) THEN
291 icomp = 1
292 jcheck = 1
293 ELSE
294 iend_finding = 1
295 ENDIF
296 ENDIF ! JCHECK
297 EXIT
298 ENDIF ! IFIT_SUCCESS
299 ENDDO ! WHILE ( ISTART < NSTART )
300 IF(iswitch == 0 .AND. ifit_success == 0 ) THEN
301 IF (jcheck == 2) THEN
302 jcheck = 1
303 iend_iter = 2
304 ELSEIF(jcheck == 1) THEN
305 IF(lawid == 2) THEN
306 lawid = 1
307 jcheck = 2
308 iend_iter =1
309 ELSEIF(nmual < 10) THEN
310 nmual = min(10, nmual + 2)
311 jcheck = 2
312 iend_iter = 1
313 ELSEIF(icomp == 0) THEN
314 icomp = 1
315 jcheck = 1
316 ELSE
317 iend_finding = 1
318 ENDIF
319 ENDIF ! JCHECK
320 ENDIF
321 IF(ifit_success == 1 .OR. iend_finding == 1) EXIT
322 ENDDO ! IFIT_SUCCESS
323 DEALLOCATE (sig)
324
325 IF (ifit_success == 0) THEN
326 IF (npsaved == 0) THEN
327 CALL ancmsg(msgid=901,
328 . msgtype=msgerror,
329 . anmode=aninfo,
330 . i1=id,
331 . c1=titr)
332 ENDIF
333 ENDIF
334C
335 nmual = nmuals
336 DO i=1,10
337 mual(i)=as(i)
338 ENDDO
339
340 IF (lawid == 2) THEN
341 DO i=5,10
342 mual(i) = zero
343 ENDDO
344 ENDIF
345 IF(is_encrypted == 0)THEN
346 WRITE(iout,'(//6X,A,/)')'FITTING RESULT COMPARISON:'
347 WRITE(iout,'(6X,A,/)')'uniaxial test data'
348 WRITE(IOUT,'(a20,5x,a20,a30)') 'nominal strain',
349 * 'nominal stress(test)', 'nominal stress(radioss)'
350 DO I=1,NPT
351 CALL OGDEN0(STRETCH(I), AS, YOGD, NMUALS)
352 WRITE(IOUT, '(1f18.6, 3x,1f20.13, 6x, 1f20.13)')
353 * STRETCH(I)-ONE,Y(I),YOGD
354 ENDDO
355 ENDIF
356!
357 WRITE(IOUT, '(a)') '-------------------------------------------'
358 WRITE(IOUT, '(a,f10.2,a)') 'averaged error of fitting : ',
359 . ERRAVE_MIN*100.0, '%'
360C-----------
361 RETURN
362 END
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
law69_nlsqf_auto
subroutine law69_nlsqf_auto(lawid, stretch, y, ma, nmual, npt, mual, icheck, nstart, errtol, id, titr, is_encrypted)
Definition law69_nlsqf_auto.F:39
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
mrqmin
subroutine mrqmin(x, y, sig, ndata, a, ia, ma, covar, alpha, nca, errnow, ifuncs, gamma, iret, icomp, mmax, atry)
Definition nlsqf.F:45
law69_check
subroutine law69_check(lawid, a, nmual, icheck, mu_incr, ivalid)
@purpose check if the guess of material properties is valid for LM optimization
Definition nlsqf.F:1180
law69_guess_bounds
subroutine law69_guess_bounds(lawid, icheck, x, y, npt, nmual, mcof_min, mcof_max, icomp)
@purpose Esstimate lwer/upper bounds of mu_i and alpha_i which will be used to generate initial guess
Definition nlsqf.F:1275
law69_guess1
subroutine law69_guess1(lawid, nmual, icheck, mu_incr, irnd1, a0, ia, mcof_min, mcof_max, nguess)
@purpose get one guess for starting Levenberg-Marquardt method
Definition nlsqf.F:1111
ogden0
subroutine ogden0(x, a, y, na)
@purpose calculate stress from stretch for OGDEN (uniaxial tension)
Definition nlsqf.F:473
radioss
program radioss
Definition radioss.F:34
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