OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
nlsqf.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/.
23C IDEBUG=1 -> output debug information into .out file
24#define IDEBUG 0
25C IOCSV=1 -> output fitting to curv.csv
26#define IOCSV 0
27!> @param[in] IFUNCS SUBROUTINE handle/ID for evaluation : 1=OGDEN0/OGDEN1
28!> @param[out] A updated coeeficients
29!> @param[out] IRET return code
30!> 0 : optimization succeeded;
31!> 1 : optimization failed;
32!||====================================================================
33!|| mrqmin ../starter/source/materials/tools/nlsqf.F
34!||--- called by ------------------------------------------------------
35!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
36!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
37!||--- calls -----------------------------------------------------
38!|| covsrt ../starter/source/materials/tools/nlsqf.F
39!|| inversion ../starter/source/materials/tools/nlsqf.F
40!|| mrqcof ../starter/source/materials/tools/nlsqf.F
41!||====================================================================
42 SUBROUTINE mrqmin(X,Y,SIG,NDATA,A,IA,
43 . MA,COVAR,ALPHA,NCA, ERRNOW,
44 . IFUNCS,GAMMA,IRET,ICOMP,MMAX,ATRY )
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "units_c.inc"
53C-----------------------------------------------
54C D u m m y A r g u m e n t s
55C-----------------------------------------------
56 INTEGER , INTENT(IN) :: MMAX
57 my_real , INTENT(INOUT) :: ATRY(MMAX)
58C----------------------------------------------------------------
59C L O C A L V A R I B L E S
60C----------------------------------------------------------------
61 INTEGER MA,NCA,NDATA,IA(MA) ,ICOMP
62 my_real
63 . gamma,errnow,a(ma),alpha(nca,nca),
64 . covar(nca,nca),x(ndata),y(ndata),sig(ndata)
65 INTEGER J,K,L,MFIT ,NMAX
66 PARAMETER (NMAX = 20) ! NMAX = MMAX
67 my_real errpre,beta(nmax),da(nmax)
68 SAVE errpre,beta,da,mfit
69 INTEGER IFUNCS, IRET
70C=======================================================================
71
72 iret = 0
73
74 IF (gamma < zero) THEN
75 mfit=0
76 DO j=1,ma
77 IF (ia(j)/=0) mfit=mfit+1
78 ENDDO
79! in Numerical Recript, 0.001 is used, however from my tests, the
80! difference is small (0.01 vs 0.001)
81 gamma=0.01
82 CALL mrqcof(x,y,sig,ndata,
83 . a,ia,ma,alpha,beta,nca,errnow,
84 . ifuncs,icomp)
85
86 errpre=errnow
87 DO j=1,ma
88 atry(j)=a(j)
89 ENDDO
90 ENDIF
91C
92 DO j=1,mfit
93 DO k=1,mfit
94 covar(j,k)=alpha(j,k)
95 ENDDO
96 covar(j,j)=alpha(j,j)*(1.+gamma)
97 da(j)=beta(j)
98 ENDDO
99C
100 CALL inversion(covar,mfit,nca,da,1,1,iret)
101 IF (iret /= 0) THEN
102 IF (idebug > 0) THEN
103 WRITE(*,*) ' IRET = ', iret
104 DO j=1,ma
105 WRITE(iout, *) 'J, A(J) = ', j, a(j)
106 ENDDO
107 ENDIF
108 RETURN
109 ENDIF
110
111 IF (gamma == zero) THEN
112 CALL covsrt(covar,nca,ma,ia,mfit)
113 RETURN
114 ENDIF
115C
116 j=0
117 DO l=1,ma
118 IF(ia(l)/=0) THEN
119 j=j+1
120 atry(l)=a(l)+da(j)
121 ENDIF
122 ENDDO
123
124 IF (idebug > 0) THEN
125 WRITE(iout,*) __file__,__line__
126 DO j=1, ma
127 write(iout,*) 'J,ATRY(J) = ', j, atry(j)
128 ENDDO
129 ENDIF
130
131c limit ATRY (ALPHA) to avoid too big power (overflow)
132c |ALPHA_i| <= 20.0
133 DO j=1, ma-1, 2
134 if (atry(j+1) > twenty) then
135 IF (idebug > 0) THEN
136 write(iout,*) 'Setting Big ALPHA to ',twenty
137 write(iout,*) ' J, ATRY(J) = ',j+1,atry(j+1)
138 ENDIF
139 atry(j+1) = twenty
140 ELSEIF (atry(j+1) < -twenty) THEN
141 IF (idebug > 0) THEN
142 write(iout,*) 'Setting Big ALPHA to -',twenty
143 write(iout,*) ' J, ATRY(J) = ',j+1,atry(j+1)
144 ENDIF
145 atry(j+1) = -twenty
146 ENDIF
147 ENDDO
148
149C
150 CALL mrqcof(x,y,sig,ndata,
151 . atry,ia,ma,covar,da,nca,errnow,
152 . ifuncs,icomp)
153C
154 IF (errnow < errpre) THEN
155 gamma=0.1*gamma
156 errpre=errnow
157 DO j=1,mfit
158 DO k=1,mfit
159 alpha(j,k)=covar(j,k)
160 ENDDO
161 beta(j)=da(j)
162 ENDDO
163 DO l=1,ma
164 a(l)=atry(l)
165 ENDDO
166 ELSE
167 gamma=10.*gamma
168 errnow=errpre
169 ENDIF
170C-----------
171 RETURN
172 END SUBROUTINE mrqmin
173C-----------
174C
175C
176!||====================================================================
177!|| mrqcof ../starter/source/materials/tools/nlsqf.F
178!||--- called by ------------------------------------------------------
179!|| mrqmin ../starter/source/materials/tools/nlsqf.F
180!||--- calls -----------------------------------------------------
181!|| my_exit ../starter/source/output/analyse/analyse.c
182!|| ogden0 ../starter/source/materials/tools/nlsqf.F
183!|| ogden1 ../starter/source/materials/tools/nlsqf.F
184!||====================================================================
185 SUBROUTINE mrqcof(X,Y,SIG, NDATA,
186 * A,IA,MA,ALPHA,BETA,NALP,
187 * ERRNOW,IFUNCS,ICOMP)
188C-----------------------------------------------
189C I m p l i c i t T y p e s
190C-----------------------------------------------
191#include "implicit_f.inc"
192C-----------------------------------------------
193C D u m m y A r g u m e n t s
194C-----------------------------------------------
195 INTEGER MA,NALP,NDATA,IA(MA),ICOMP
196 my_real
197 . ERRNOW,A(MA),ALPHA(NALP,NALP),BETA(MA),
198 . X(NDATA),Y(NDATA),SIG(NDATA)
199 INTEGER IFUNCS
200C----------------------------------------------------------------
201C L O C A L V A R I B L E S
202C----------------------------------------------------------------
203 INTEGER MFIT,I,J,K,L,M,MMAX
204 PARAMETER (MMAX=20)
205 my_real
206 . dy,wt,ymod,dyda(mmax)
207
208 my_real y_min
209
210 y_min = em3
211
212C=======================================================================
213 mfit=0
214 DO 11 j=1,ma
215 IF (ia(j)/=0) mfit=mfit+1
21611 CONTINUE
217 DO 13 j=1,mfit
218 DO 12 k=1,j
219 alpha(j,k)=0.
22012 CONTINUE
221 beta(j)=0.
22213 CONTINUE
223 errnow=zero
224 IF(icomp == 0) THEN
225 DO 16 i=1,ndata
226 IF (ifuncs == 1) THEN
227 CALL ogden0(x(i),a,ymod, ma)
228 CALL ogden1(x(i),a,dyda, ia, ma)
229 ELSE
230! TBD, UNKOWN IFUNCS, ERROR OUT!!!
231 WRITE(*,*) 'ERROR: UNKOWN IFUNCS = ', ifuncs
232 CALL my_exit(2)
233 ENDIF
234
235 dy=y(i)-ymod
236 errnow=errnow+dy*dy/(sig(i)*sig(i))
237c IF (Y(I)>EM5) DY=Y(I)-YMOD
238c IF (Y(I)<EM5) DY=EM5
239
240 j=0
241 DO 15 l=1,ma
242 IF(ia(l)/=0) THEN
243 j=j+1
244 wt=dyda(l)
245 wt=dyda(l)/(sig(i)*sig(i))
246 k=0
247 DO 14 m=1,l
248 IF(ia(m)/=0) THEN
249 k=k+1
250 alpha(j,k)=alpha(j,k)+wt*dyda(m)
251 ENDIF
25214 CONTINUE
253 beta(j)=beta(j)+dy*wt
254 ENDIF
25515 CONTINUE
25616 CONTINUE
257 ELSE
258 DO i=1,ndata
259 IF (ifuncs == 1) THEN
260 CALL ogden0(x(i),a,ymod, ma)
261 CALL ogden1(x(i),a,dyda, ia, ma)
262 ELSE
263! TBD, UNKOWN IFUNCS, ERROR OUT!!!
264 WRITE(*,*) 'ERROR: UNKOWN IFUNCS = ', ifuncs
265 CALL my_exit(2)
266 ENDIF
267
268 dy=y(i)-ymod
269 errnow = errnow + dy*dy
270 j=0
271 DO l=1,ma
272 IF(ia(l)/=0) THEN
273 j=j+1
274 wt=dyda(l)
275 k=0
276 DO m=1,l
277 IF(ia(m)/=0) THEN
278 k=k+1
279 alpha(j,k)=alpha(j,k)+wt*dyda(m)
280 ENDIF
281 ENDDO ! M
282 beta(j)=beta(j)+dy*wt
283 ENDIF
284 ENDDO ! L
285 ENDDO ! I
286 ENDIF
287 DO 18 j=2,mfit
288 DO 17 k=1,j-1
289 alpha(k,j)=alpha(j,k)
29017 CONTINUE
291
29218 CONTINUE
293 RETURN
294 END SUBROUTINE mrqcof
295C
296
297!||====================================================================
298!|| inversion ../starter/source/materials/tools/nlsqf.F
299!||--- called by ------------------------------------------------------
300!|| mrqmin ../starter/source/materials/tools/nlsqf.F
301!|| mrqmin_law92 ../starter/source/materials/mat/mat092/law92_nlsqf.F90
302!||====================================================================
303 SUBROUTINE inversion(A,N,NP,B,M,MP,IRET)
304C-----------------------------------------------
305C I M P L I C I T T Y P E S
306C-----------------------------------------------
307#include "implicit_f.inc"
308C-----------------------------------------------
309C D u m m y A r g u m e n t s
310C-----------------------------------------------
311 INTEGER M,MP,N,NP
312 my_real A(NP,NP),B(NP,MP)
313C----------------------------------------------------------------
314C L O C A L V A R I B L E S
315C----------------------------------------------------------------
316 INTEGER I,NMAX,ICOL,IROW,J,K,L,LL
317 PARAMETER (NMAX=50)
318 integer indxc(nmax),indxr(nmax),ipiv(nmax)
319 my_real
320 . big,dum,pivinv
321
322 INTEGER IRET
323
324 iret = 0
325
326 icol = -1
327 irow = -1
328
329C=======================================================================
330 DO 11 j=1,n
331 ipiv(j)=0
33211 CONTINUE
333 DO 22 i=1,n
334 big=zero
335 DO 13 j=1,n
336 IF(ipiv(j)/=1)THEN
337 DO 12 k=1,n
338 IF (ipiv(k)==0) THEN
339 IF (abs(a(j,k))>=big)THEN
340 big=abs(a(j,k))
341 irow=j
342 icol=k
343 ENDIF
344 ELSE IF (ipiv(k)>1) THEN
345 iret = __line__
346 WRITE(*,*) 'SINGULAR MATRIX IN INVERSION'
347 RETURN
348 ENDIF
34912 CONTINUE
350 ENDIF
35113 CONTINUE
352
353 IF (icol == -1) THEN
354 iret = __line__
355 WRITE(*,*) 'FATAL ERROR IN INVERSION'
356 RETURN
357 ENDIF
358
359 ipiv(icol)=ipiv(icol)+1
360 IF (irow/=icol) THEN
361 DO 14 l=1,n
362 dum=a(irow,l)
363 a(irow,l)=a(icol,l)
364 a(icol,l)=dum
36514 CONTINUE
366 DO 15 l=1,m
367 dum=b(irow,l)
368 b(irow,l)=b(icol,l)
369 b(icol,l)=dum
37015 CONTINUE
371 ENDIF
372 indxr(i)=irow
373 indxc(i)=icol
374 IF (a(icol,icol)==0.) THEN
375 iret = __line__
376 WRITE(*,*) 'SINGULAR MATRIX IN INVERSION'
377 RETURN
378 ENDIF
379 pivinv=1./a(icol,icol)
380 a(icol,icol)=1.
381 DO 16 l=1,n
382 a(icol,l)=a(icol,l)*pivinv
38316 CONTINUE
384 DO 17 l=1,m
385 b(icol,l)=b(icol,l)*pivinv
38617 CONTINUE
387 DO 21 ll=1,n
388 IF(ll/=icol)THEN
389 dum=a(ll,icol)
390 a(ll,icol)=0.
391 DO 18 l=1,n
392 a(ll,l)=a(ll,l)-a(icol,l)*dum
39318 CONTINUE
394 DO 19 l=1,m
395 b(ll,l)=b(ll,l)-b(icol,l)*dum
39619 CONTINUE
397 ENDIF
39821 CONTINUE
39922 CONTINUE
400 DO 24 l=n,1,-1
401 IF(indxr(l)/=indxc(l))THEN
402 DO 23 k=1,n
403 dum=a(k,indxr(l))
404 a(k,indxr(l))=a(k,indxc(l))
405 a(k,indxc(l))=dum
40623 CONTINUE
407 ENDIF
40824 CONTINUE
409C-----------
410 RETURN
411 END SUBROUTINE inversion
412C
413C
414!||====================================================================
415!|| covsrt ../starter/source/materials/tools/nlsqf.F
416!||--- called by ------------------------------------------------------
417!|| mrqmin ../starter/source/materials/tools/nlsqf.F
418!||====================================================================
419 SUBROUTINE covsrt(COVAR,SIZCOVAR,MA,IA,MFIT)
420C-----------------------------------------------
421C I M P L I C I T T Y P E S
422C-----------------------------------------------
423#include "implicit_f.inc"
424C----------------------------------------------------------------
425C D u m m y A r g u m e n t s
426C----------------------------------------------------------------
427 INTEGER MA,MFIT,SIZCOVAR,IA(MA)
428 my_real COVAR(SIZCOVAR,SIZCOVAR),SWAP
429C----------------------------------------------------------------
430C L O C A L V A R I B L E S
431C----------------------------------------------------------------
432 INTEGER I,J,K
433C=======================================================================
434 DO 12 I=mfit+1,ma
435 DO 11 j=1,i
436 covar(i,j)=0.
437 covar(j,i)=0.
43811 CONTINUE
43912 CONTINUE
440 k=mfit
441 DO 15 j=ma,1,-1
442 IF(ia(j)/=0)THEN
443 DO 13 i=1,ma
444 swap=covar(i,k)
445 covar(i,k)=covar(i,j)
446 covar(i,j)=swap
44713 CONTINUE
448 DO 14 i=1,ma
449 swap=covar(k,i)
450 covar(k,i)=covar(j,i)
451 covar(j,i)=swap
45214 CONTINUE
453 k=k-1
454 ENDIF
45515 CONTINUE
456C-----------
457 RETURN
458 END SUBROUTINE covsrt
459
460!> @purpose
461!> calculate stress from stretch for OGDEN (uniaxial tension)
462!> @param[in] X stretch
463!> @param[in] A(NA) material parameters [mu_1, alpha_1,mu_2, alpha_2,...]
464!> @param[out] Y stress
465!||====================================================================
466!|| ogden0 ../starter/source/materials/tools/nlsqf.F
467!||--- called by ------------------------------------------------------
468!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
469!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
470!|| mrqcof ../starter/source/materials/tools/nlsqf.F
471!||====================================================================
472 SUBROUTINE ogden0(X,A,Y,NA)
473C-----------------------------------------------
474C I m p l i c i t T y p e s
475C-----------------------------------------------
476#include "implicit_f.inc"
477C-----------------------------------------------
478C D u m m y A r g u m e n t s
479C-----------------------------------------------
480 INTEGER NA
481 my_real X, A(NA), Y
482C----------------------------------------------------------------
483C L O C A L V A R I B L E S
484C----------------------------------------------------------------
485 INTEGER I
486C=======================================================================
487 Y=zero
488 DO i=1,na-1,2
489 y = y + a(i)*(x**(-1.+a(i+1))-x**(-0.5*a(i+1)-1.))
490 ENDDO
491C-----------
492 RETURN
493 END SUBROUTINE ogden0
494
495!> @purpose
496!> calculate d(stress)/d(A) OGDEN (uniaxial tension)
497!> @param[in] X stretch
498!> @param[in] A(NA) material parameters [mu_1, alpha_1,mu_2, alpha_2,...]
499!> @param[in] IA(NA) flags if derivative needs to be calculated
500!> 1=yes, 0=no
501!> @param[out] DYDA(NA) derivatives d(stress)/d(A)
502!||====================================================================
503!|| ogden1 ../starter/source/materials/tools/nlsqf.F
504!||--- called by ------------------------------------------------------
505!|| mrqcof ../starter/source/materials/tools/nlsqf.F
506!||====================================================================
507 SUBROUTINE ogden1(X,A,DYDA,IA,NA)
508C-----------------------------------------------
509C I m p l i c i t T y p e s
510C-----------------------------------------------
511#include "implicit_f.inc"
512C-----------------------------------------------
513C D u m m y A r g u m e n t s
514C-----------------------------------------------
515 INTEGER NA
516 INTEGER IA(NA)
517 my_real X, A(NA),DYDA(NA)
518 my_real TMP1, TMP2
519C----------------------------------------------------------------
520C L O C A L V A R I B L E S
521C----------------------------------------------------------------
522 INTEGER I
523C=======================================================================
524 DO I=1,na-1,2
525 tmp1 = x**(-1.+a(i+1))
526 tmp2 = x**(-0.5*a(i+1)-1.)
527 IF (ia(i) /= 0) THEN
528 dyda(i) = tmp1 - tmp2
529 ENDIF
530
531 IF (ia(i+1) /= 0) THEN
532 dyda(i+1) = a(i)*log(x)*(tmp1 + 0.5*tmp2)
533 ENDIF
534 ENDDO
535
536 RETURN
537 END SUBROUTINE ogden1
538
539!> @param[in] LAWID law_ID on /MAT/LAW69 (1=Ogden, 2=Mooney-Rivlin)
540!> @param[in] ICHECK checking level
541!> <=0 no constraint (IVALID always return as 1)
542!> 1 SUM( mu(i) * alpha(i) ) > 0.0
543!> 2 mu(i) * alpha(i) > 0.0
544!> 3 Try ICHECK=2 at first, if it fails in fitting, switch to ICHECK=1, and try again
545!> @param[in] ERRTOL If ERRAVE < ERRTOL, data fitting converges.
546!> ERRAVE = ( SUM [ ABS ( ( Y_inp-Y_fit) / Y_inp ) ) / NPT
547!> @param[out] MUAL(1:NMUAL) optimized material properties
548!||====================================================================
549!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
550!||--- called by ------------------------------------------------------
551!|| law69_upd ../starter/source/materials/mat/mat069/law69_upd.F
552!||--- calls -----------------------------------------------------
553!|| ancmsg ../starter/source/output/message/message.F
554!|| law69_check ../starter/source/materials/tools/nlsqf.f
555!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
556!|| law69_guess_bounds ../starter/source/materials/tools/nlsqf.F
557!|| mrqmin ../starter/source/materials/tools/nlsqf.F
558!|| my_exit ../starter/source/output/analyse/analyse.c
559!|| ogden0 ../starter/source/materials/tools/nlsqf.F
560!||--- uses -----------------------------------------------------
561!|| message_mod ../starter/share/message_module/message_mod.F
562!||====================================================================
563 SUBROUTINE law69_nlsqf(LAWID,STRETCH,Y,MA,NMUAL,NPT,MUAL,
564 $ ICHECK, NSTART, ERRTOL,ID,TITR,IS_ENCRYPTED)
565 USE message_mod
566 USE names_and_titles_mod , ONLY : nchartitle
567C-----------------------------------------------
568C I m p l i c i t T y p e s
569C-----------------------------------------------
570#include "implicit_f.inc"
571C-----------------------------------------------
572C C o m m o n B l o c k s
573C-----------------------------------------------
574#include "units_c.inc"
575C-----------------------------------------------
576 INTEGER MAXA
577 PARAMETER (MAXA=20)
578
579 integer lawid, npt, ma, icheck, nstart
580 INTEGER , INTENT(IN) :: IS_ENCRYPTED
581 my_real errtol
582 INTEGER I,NONZERO(MAXA),IDUM,ITER,ICOUNT,J,K,NMUAL
583 my_real GAMMA,ERRNOW,GASDEV,ERRPRE,
584 . stretch(npt),mual(10)
585 my_real a(maxa),covar(maxa,maxa),alpha(maxa,maxa),y(npt)
586 my_real mcof_min(maxa), mcof_max(maxa)
587 my_real yogd,xogd
588
589 INTEGER ID
590 CHARACTER(LEN=NCHARTITLE) :: TITR
591 INTEGER MINITER_LM, MAXITER_LM
592 DATA miniter_lm /5/
593 SAVE miniter_lm
594
595 DATA maxiter_lm /100/
596 SAVE maxiter_lm
597
598 ! if ( (ERRNOW-ERRPRE) / ERRPRE ) < EPS_LM, CNT_HIT_EPS_LM = CNT_HIT_EPS_LM + 1,
599 ! if CNT_HIT_EPS_LM == LMT_HIT_EPS_LM, it converges
600 my_real eps_lm
601 DATA eps_lm /1e-3/
602 SAVE eps_lm
603
604 INTEGER CNT_HIT_EPS_LM
605 INTEGER LMT_HIT_EPS_LM
606 DATA lmt_hit_eps_lm /2/
607 SAVE lmt_hit_eps_lm
608
609 INTEGER LMSTOP
610
611 INTEGER IFUNCS
612
613 INTEGER NGUESS
614 my_real a0(maxa)
615 my_real errmin, errmax, errave, errave_min, err
616
617 INTEGER ISTART, NPSAVED, IVALID,ICOMP, MMAX
618 parameter(mmax =20)
619 my_real atry(mmax)
620
621 INTEGER ICURV
622 LOGICAL lopened
623
624c during LM optimization, if the objective is not improved,
625c GAMMA will be increased. We should terminate the iteration once
626c GAMMA becomes very huge.
627 my_real gamma_stop
628 data gamma_stop /1e15/
629 save gamma_stop
630
631c if check the validity of initial guess (0=no, 1=yes)
632c we don't want to check, because invalid initial point could converge to valid point.
633 INTEGER ICHECK_GUESS
634 data icheck_guess /0/
635 save icheck_guess
636
637 INTEGER JCHECK, IFIT_SUCCESS
638
639c if enforce mu(i) < mu(i+1) in generating initial guess
640c we don't need this anymore since we are using random numbers instead of loop through all
641c the combinations
642 INTEGER MU_INCR_GUESS
643 data MU_INCR_GUESS /0/
644 save MU_INCR_GUESS
645
646 my_real MAX_ABS_YI
647 my_real small_fac_abs_yi, small_abs_yi
648
649 my_real, DIMENSION(:), ALLOCATABLE :: sig
650 my_real spready
651 INTEGER IRET
652 INTEGER IRND1
653! ------------------------------------------------------------------------------
654 ALLOCATE (sig(1:npt))
655!
656 !
657 IF ( maxa < ma ) THEN
658 WRITE(*,*) 'ERROR, MAXA < MA'
659 WRITE(*,*) __file__,__line__
660 CALL my_exit(2)
661 ENDIF
662 icomp=0
663 IF(icheck < 0) icomp= 1
664 icheck = abs(icheck)
665 ! IF ABS(Y(I)) < SMALL_ABS_YI, use SMALL_ABS_YI to avoid
666 ! unnecessary numerical issues when avoid divided by small value
667
668 max_abs_yi = zero
669 DO i = 1, npt
670 max_abs_yi = max( max_abs_yi, abs(y(i)) )
671 ENDDO
672
673 small_fac_abs_yi = em3
674
675 small_abs_yi = max_abs_yi * small_fac_abs_yi
676
677 IF (idebug > 0) THEN
678 WRITE(iout, *) ' MAX_ABS_YI, SMALL_FAC_ABS_YI, SMALL_ABS_YI = ',
679 $ max_abs_yi, small_fac_abs_yi, small_abs_yi
680 ENDIF
681
682 spready=0.01
683 IF(icomp == 0) THEN
684 DO j=1,npt
685 sig(j)=spready*max(small_abs_yi, abs(y(j)) )
686 IF (idebug > 0) THEN
687 WRITE(iout, *) 'J, SIG(J) = ', j, sig(j)
688 ENDIF
689 ENDDO
690 ELSE
691 DO j=1,npt
692 sig(j) = y(j)
693 IF(sig(j) == zero) sig(j) = em15
694 IF (idebug > 0) THEN
695 WRITE(iout, *) 'J, SIG(J) = ', j, sig(j)
696 ENDIF
697 ENDDO
698 ENDIF
699
700 IF (maxa < ma) THEN
701! TBD_KANG
702 WRITE(*,*) 'ERROR: MAXA < MA'
703 WRITE(*,*) ' MAXA = ', maxa
704 WRITE(*,*) ' MA = ', maxa
705 WRITE(*,*) ' FILE = ', __file__
706 WRITE(*,*) ' LINE = ', __line__
707 CALL my_exit(2)
708 ENDIF
709
710C=======================================================================
711 DO i=1,nmual
712 nonzero(i)=1
713 ENDDO
714
715 IF (lawid == 1) then ! OGDEN
716 DO i=nmual+1,10
717 nonzero(i)=0
718 a(i)=zero
719 ENDDO
720
721 ifuncs = 1 ! SUBROUTINE OGDEN0/OGDEN1
722
723 ELSEIF (lawid == 2) then ! Mooney-Rivlin
724 DO i=nmual+1,ma
725 nonzero(i)=0
726 ENDDO
727 nonzero(2)=0
728 nonzero(4)=0
729
730 DO i=nmual+1,ma
731 a(i)=zero
732 ENDDO
733
734 ifuncs = 1 ! SUBROUTINE OGDEN0/OGDEN1
735
736 ENDIF
737
738 jcheck =icheck
739 IF (jcheck == 3) THEN
740 jcheck = 2
741 ENDIF
742
74399 CONTINUE ! location to start optimization with different JCHECK
744
745 ifit_success = 0
746
747 istart = 0
748 npsaved = 0
749
750c calculate MCOF_MIN, MCOF_MAX
751 CALL law69_guess_bounds(lawid, jcheck, stretch, y, npt,
752 $ nmual, mcof_min, mcof_max,icomp)
753
754c initialize random number generator in LAW69_GUESS1
755 irnd1 = 205
756 atry(1:mmax) = zero ! Initialization of ATRY
757 DO 111 WHILE ( istart < nstart )
758c get one guess point A0(1:NMUAL)
759 CALL law69_guess1(
760 $ lawid, nmual, icheck_guess, mu_incr_guess,irnd1,
761 $ a0, nonzero, mcof_min, mcof_max, nguess)
762
763 IF (nguess == 0) THEN
764 GOTO 112
765 ENDIF
766
767C calculate averaged ERROR before LM
768 errave = 0.0
769 IF(icomp == 0) THEN
770 DO i=1,npt
771 CALL ogden0(stretch(i), a0, yogd, nmual)
772 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
773 errave = errave + err
774 ENDDO
775 ELSE
776 DO i=1,npt
777 CALL ogden0(stretch(i), a0, yogd, nmual)
778 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
779 errave = errave + err
780 ENDDO
781 ENDIF
782C
783 errave = errave / (1.0 * npt)
784
785 istart = istart + 1
786
787 IF (idebug > 0) THEN
788 WRITE(iout,*) ' ISTART = ', istart
789 WRITE(iout,*) ' Before LM optimization ...'
790 DO i = 1, nmual
791 WRITE(iout, *) ' I, A0(I) = ', i, a0(i)
792 ENDDO
793 WRITE(iout,*) ' LM0, ERRAVE = ', errave
794 ENDIF
795
796! start loop of LM optimization
797 iter = 0
798 cnt_hit_eps_lm = 0
799
800 iter = iter + 1
801 gamma=-1
802 CALL mrqmin(stretch,y,sig,npt,a0,nonzero,
803 $ nmual,covar,alpha,ma,errnow,
804 $ ifuncs,gamma,iret,icomp,mmax,atry)
805 IF (iret /= 0) GOTO 111
806
807 IF (idebug > 0) THEN
808 errave = 0.0
809 IF(icomp == 0) THEN
810 DO i=1,npt
811 CALL ogden0(stretch(i), a0, yogd, nmual)
812 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
813 errave = errave + err
814 ENDDO
815 ELSE
816 DO i=1,npt
817 CALL ogden0(stretch(i), a0, yogd, nmual)
818 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
819 errave = errave + err
820 ENDDO
821 ENDIF
822 errave = errave / (1.0 * npt)
823
824 WRITE(iout, '(A,I4, 3E16.8)')
825 $ 'ITER, ERRNOW, GAMMA, ERRAVE = ',
826 $ iter, errnow, gamma, errave
827 DO i = 1, nmual
828 WRITE(iout, *) ' - I, A0(I) = ', i, a0(i)
829 ENDDO
830 ENDIF
831
83221 CONTINUE
833 errpre=errnow
834
835 iter = iter + 1
836 CALL mrqmin(stretch,y,sig,npt,a0,nonzero,
837 $ ma,covar,alpha, ma, errnow,
838 $ ifuncs,gamma,iret,icomp,mmax,atry)
839 IF (iret /= 0) GOTO 111
840
841 IF (idebug == 1) THEN
842 errave = zero
843 IF(icomp == 0) THEN
844 DO i=1,npt
845 CALL ogden0(stretch(i), a0, yogd, nmual)
846 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
847 errave = errave + err
848 ENDDO
849 ELSE
850 DO i=1,npt
851 CALL ogden0(stretch(i), a0, yogd, nmual)
852 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
853 errave = errave + err
854 ENDDO
855 ENDIF
856 errave = errave / (1.0 * npt)
857
858 WRITE(iout, '(A,I4, 5E16.8)')
859 $ 'ITER, ERRNOW, GAMMA, ERRAVE, ERRNOW-ERRPRE,'//
860 $ '(ERRNOW-ERRPRE)/ERRPRE = ',
861 $ iter, errnow, gamma, errave, errnow-errpre,
862 $ (errnow-errpre)/errpre
863 DO i = 1, nmual
864 WRITE(iout, *) ' - I, A0(I) = ', i, a0(i)
865 ENDDO
866 ENDIF
867
868c IF A0(J) is too small, restart from next initial guess
869 DO j = 1, nmual
870 IF ( abs( a0(j) ) < 1e-20 ) THEN
871 goto 111 ! restart from next initial guess
872c WRITE(IOUT, *) ' A is too small, and could result in error'
873c WRITE(IOUT,*) 'J, A0(J) = ', J, A0(J)
874c STOP
875 ENDIF
876 ENDDO
877
878c check convergence of LM optimization
879 lmstop = 0
880 IF (idebug > 0) THEN
881 WRITE(iout, *) ' ERRNOW/(1.0*NPT) = ', errnow/(1.0*npt)
882 WRITE(iout, *) ' ERRNOW < ERRPRE = ', (errnow < errpre)
883 ENDIF
884
885 IF ( iter > miniter_lm ) THEN
886 IF (errnow < errpre) THEN
887 IF ( abs(errpre) > zero) THEN
888 IF ( abs( (errnow-errpre)/ errpre ) < eps_lm) THEN
889 cnt_hit_eps_lm = cnt_hit_eps_lm + 1
890
891 IF (idebug > 0) THEN
892 WRITE(iout,*)
893 $ ' CNT_HIT_EPS_LM, ABS((ERRNOW-ERRPRE)/ERRPRE) = ',
894 $ cnt_hit_eps_lm, abs( (errnow-errpre)/ errpre )
895 ENDIF
896
897 IF ( cnt_hit_eps_lm >= lmt_hit_eps_lm ) THEN
898 IF (idebug > 0) THEN
899 WRITE(iout,*) 'STOP AT ', __line__
900 ENDIF
901 lmstop = 1
902 ENDIF
903 ENDIF
904 ENDIF
905 ELSEIF (iter >= maxiter_lm .OR. gamma >= gamma_stop) THEN
906 IF (idebug > 0) THEN
907 WRITE(iout,*) 'STOP AT ', __line__
908 ENDIF
909 lmstop = 1
910 ENDIF
911 ENDIF
912
913 IF (lmstop== 0) THEN
914 GOTO 21
915 ENDIF
916! end loop of LM optimization
917
918 errave = 0.0
919 IF(icomp == 0) THEN
920 DO i=1,npt
921 CALL ogden0(stretch(i), a0, yogd, nmual)
922 err = abs(y(i) - yogd) / max(small_abs_yi, abs(y(i)))
923 errave = errave + err
924 ENDDO
925 ELSE
926 DO i=1,npt
927 CALL ogden0(stretch(i), a0, yogd, nmual)
928 err = abs(y(i) - yogd) / max(em15,abs(y(i)))
929 errave = errave + err
930 ENDDO
931 ENDIF
932 errave = errave / (1.0 * npt)
933
934 IF (idebug > 0) THEN
935 WRITE(iout,*) ' After LM optimization ...'
936 DO i = 1, ma
937 WRITE(iout, *) ' I, A0(I) = ', i, a0(i)
938 ENDDO
939 WRITE(iout,*) ' LM1, ERRAVE = ', errave
940 ENDIF
941
942 CALL law69_check(lawid, a0, nmual, jcheck, 0, ivalid)
943 IF (ivalid > 0) THEN
944 IF ( npsaved==0 .OR.
945 $ ( npsaved>0 .AND. errave<errave_min) ) THEN
946 npsaved = npsaved + 1
947 errave_min = errave
948 DO i = 1, nmual
949 a(i) = a0(i)
950 ENDDO
951 ENDIF
952 ELSE
953 IF (idebug > 0) THEN
954 WRITE(iout,*) __file__,__line__
955 WRITE(iout,*) ' LM converges to invalid point'
956 ENDIF
957 ENDIF
958
959 IF (npsaved > 0) THEN
960 IF (idebug > 0) THEN
961 WRITE(*, *) ' ISTART, NPSAVED, ERRAVE_MIN = ',
962 $ istart, npsaved, errave_min
963 WRITE(iout, *) ' ISTART, NPSAVED, ERRAVE_MIN = ',
964 $ istart, npsaved, errave_min
965 ENDIF
966
967 IF ( errave_min < errtol ) THEN
968 ifit_success = 1
969 GOTO 112
970 ENDIF
971 ELSE
972 IF (idebug > 0) THEN
973 WRITE(*, *) ' ISTART, NPSAVED ',
974 $ istart, npsaved
975 WRITE(iout, *) ' ISTART, NPSAVED ',
976 $ istart, npsaved
977 ENDIF
978 ENDIF
979
980111 CONTINUE ! WHILE ( ISTART < NSTART )
981
982112 CONTINUE
983
984 IF (ifit_success == 0 .AND. icheck == 3) THEN
985 IF (jcheck == 2) THEN
986 jcheck = 1
987 IF (idebug > 0) THEN
988 IF (npsaved > 0) THEN
989 WRITE(*,'(A)')
990 $ '*** Curve fitting result of /MAT/LAW69 with'
991 $ //' ICHECK=2 is not satisfactory.'
992 WRITE(*,'(A,F10.2,A)') 'ERRAVE = ', errave_min*100.0,'%'
993 WRITE(*,'(A)') ' Switching to ICHECK=1 and trying again.'
994
995 ELSE
996 WRITE(*,'(A)')
997 $ '*** Failed to fit /MAT/LAW69 with ICHECK=2.'
998 WRITE(*,'(A)') ' Switching to ICHECK=1 and trying again.'
999 ENDIF
1000 ENDIF
1001 GOTO 99
1002 ENDIF
1003 ENDIF
1004
1005 DEALLOCATE (sig)
1006
1007 IF (ifit_success == 0) THEN
1008 IF (npsaved == 0) THEN
1009 CALL ancmsg(msgid=901,
1010 . msgtype=msgerror,
1011 . anmode=aninfo,
1012 . i1=id,
1013 . c1=titr)
1014 ENDIF
1015 ENDIF
1016C
1017 DO i=1,10
1018 mual(i)=a(i)
1019 ENDDO
1020
1021 IF (lawid == 2) THEN
1022 DO i=5,10
1023 mual(i) = zero
1024 ENDDO
1025 ENDIF
1026 IF(is_encrypted == 0)THEN
1027 WRITE(iout,'(//6X,A,/)')'FITTING RESULT COMPARISON:'
1028 WRITE(iout,'(6X,A,/)')'UNIAXIAL TEST DATA'
1029 WRITE(iout,'(a20,5x,a20,a30)') 'NOMINAL STRAIN',
1030 * 'NOMINAL STRESS(TEST)', 'NOMINAL STRESS(RADIOSS)'
1031 ENDIF
1032
1033c output curcves to .csv format to simplify visualization
1034 icurv = 0
1035 IF (iocsv> 0) THEN
1036 DO icurv = 25, 99
1037 inquire (unit=icurv, opened=lopened)
1038 if (.not. lopened) goto 77
1039 ENDDO
104077 CONTINUE
1041 OPEN(unit=icurv, file='curv.csv')
1042 WRITE(icurv,'(A)') 'NOMINAL STRAIN, NOMINAL STRESS(TEST), '//
1043 $ 'NOMINAL STRESS(RADIOSS)'
1044 ENDIF
1045 IF(is_encrypted == 0)THEN
1046 DO i=1,npt
1047 CALL ogden0(stretch(i), a, yogd, nmual)
1048!! WRITE(IOUT,'(F18.4,F20.4,F28.4)')
1049 WRITE(iout, '(1F18.6, 3X,1F20.13, 6X, 1F20.13)')
1050 * stretch(i)-one,y(i),yogd
1051 IF (icurv > 0) THEN
1052 WRITE(icurv, '(F18.4, A, F18.4, A, F18.4)')
1053 $ stretch(i)-one,',',y(i),',',yogd
1054 ENDIF
1055 ENDDO
1056 ENDIF
1057
1058 WRITE(iout, *) ''
1059 WRITE(iout, '(A)') '-------------------------------------------'
1060 WRITE(iout, '(A,F10.2,A)') 'AVERAGED ERROR OF FITTING : ',
1061 $ errave_min*100.0, '%'
1062
1063c output curcves to .csv format to simplify visualization
1064 IF ( icurv > 0) THEN
1065 CLOSE(icurv)
1066 ENDIF
1067C-----------
1068 RETURN
1069 END
1070
1071c if abs(value) < small, it means value is very close to zero,
1072c reset it to be small or -small to avoid potential numerical issue
1073!||====================================================================
1074!|| zero2small ../starter/source/materials/tools/nlsqf.F
1075!||--- called by ------------------------------------------------------
1076!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
1077!||====================================================================
1078 SUBROUTINE zero2small(value, small)
1079#include "implicit_f.inc"
1080 my_real value, small
1081
1082 if ( abs(value) < small ) then
1083 if (value >= 0.0) then
1084 value = small
1085 else
1086 value = -small
1087 endif
1088 endif
1089 END SUBROUTINE zero2small
1090
1091!> @purpose
1092!> get one guess for starting Levenberg-Marquardt method
1093!> @param[in] LAWID material type(1=Ogden, 2=Mooney-Rivlin)
1094!> @param[out] NGUESS Number of returned guess
1095!> 0 = no return,already go through all combinations,
1096!> 1 = return one guess point successfully
1097!> @param[out] A0(NMUAL) returned guess
1098!||====================================================================
1099!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
1100!||--- called by ------------------------------------------------------
1101!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
1102!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
1103!||--- calls -----------------------------------------------------
1104!|| law69_check ../starter/source/materials/tools/nlsqf.F
1105!|| rand1 ../starter/source/materials/tools/nlsqf.F
1106!|| zero2small ../starter/source/materials/tools/nlsqf.F
1107!||====================================================================
1108 SUBROUTINE law69_guess1(
1109 $ LAWID, NMUAL, ICHECK, MU_INCR, IRND1,
1110 $ A0, IA, MCOF_MIN, MCOF_MAX, NGUESS)
1111#include "implicit_f.inc"
1112C-----------------------------------------------
1113C A n a l y s e M o d u l e
1114C-----------------------------------------------
1115
1116 REAL RAND1
1117 EXTERNAL RAND1
1118
1119 INTEGER LAWID, NMUAL, ICHECK, MU_INCR, IRND1, NGUESS
1120 my_real A0(NMUAL)
1121 INTEGER IA(NMUAL)
1122
1123 my_real small
1124 data small /1e-15/
1125 save small
1126
1127 my_real mcof_min(nmual) ! lower bounds of material parameters
1128 my_real mcof_max(nmual) ! upper bounds of material parameters
1129
1130 real*4 randnum
1131
1132 INTEGER I, IVALID
1133
1134 nguess = 0
1135
1136 IF (lawid == 2) THEN
1137 a0(2) = 2.0
1138 a0(4) = -2.0
1139 ENDIF
1140
1141 DO WHILE (.true.)
1142 DO i = 1, nmual
1143 IF (ia(i) > 0) THEN
1144 randnum = rand1(irnd1)
1145 a0(i) = (1.0-randnum)*mcof_min(i) + randnum*mcof_max(i)
1146 CALL zero2small(a0(i), small)
1147 ENDIF
1148 ENDDO
1149
1150 CALL law69_check( lawid, a0, nmual,
1151 $ icheck, mu_incr, ivalid)
1152 IF (ivalid == 1) THEN
1153 nguess = 1
1154 GOTO 299
1155 ENDIF
1156 ENDDO
1157
1158299 CONTINUE
1159
1160 END SUBROUTINE law69_guess1
1161
1162!> @purpose
1163!> check if the guess of material properties is valid for LM optimization
1164!> @param[in] ICHECK checking level
1165!> <=0 no constraint (IVALID always return as 1)
1166!> 1 SUM( mu(i) * alpha(i) ) > 0.0
1167!> 2 mu(i) * alpha(i) > 0.0
1168!> @param[in] MU_INCR if condition mu(i) <= mu(i+1) is checked
1169!> 0=no, 1=yes
1170!> @param[out] IVALID 0=not valid, 1=valid
1171!||====================================================================
1172!|| law69_check ../starter/source/materials/tools/nlsqf.F
1173!||--- called by ------------------------------------------------------
1174!|| law69_guess1 ../starter/source/materials/tools/nlsqf.F
1175!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
1176!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
1177!||====================================================================
1178 SUBROUTINE law69_check(LAWID, A, NMUAL,
1179 $ ICHECK, MU_INCR, IVALID)
1180#include "implicit_f.inc"
1181#include "units_c.inc"
1182 INTEGER LAWID, NMUAL, ICHECK, MU_INCR, IVALID
1183 my_real A(NMUAL)
1184
1185c local vars
1186 INTEGER J
1187 my_real SUM_MU_X_AF
1188
1189 IF (IDEBUG > 0) THEN
1190 WRITE(iout,*) 'LAW69_CHECK ...'
1191 DO j = 1, nmual
1192 WRITE(iout,*) 'J, A(J) = ', j, a(j)
1193 ENDDO
1194 ENDIF
1195
1196 ivalid = 1
1197 IF (icheck == 0) GOTO 49
1198
1199 IF (mu_incr > 0) THEN
1200 IF (lawid == 1) THEN
1201 IF (nmual >= 4) THEN
1202 DO j = 1, nmual/2-1
1203 IF ( a(2*j-1) > a(2*(j+1)-1) ) THEN
1204 IF (idebug > 0) THEN
1205 WRITE(iout,*) 'J,2*J-1,2*(J+1)-1 = ',
1206 $ j,2*j-1,2*(j+1)-1
1207 WRITE(iout,*) ' A(2*J-1), A(2*(J+1)-1) = ',
1208 $ a(2*j-1), a(2*(j+1)-1)
1209 WRITE(iout,*) ' MU(i) > MU(i+1) '
1210 ENDIF
1211 ivalid = 0
1212 goto 49
1213 ENDIF
1214 ENDDO
1215 ENDIF
1216 ENDIF
1217 ENDIF
1218
1219 IF (icheck == 1) THEN
1220C CHECK IF SUM(MU_I * ALPHA_I) > 0.0
1221 sum_mu_x_af = 0.0
1222 DO j = 1, nmual/2
1223 sum_mu_x_af = sum_mu_x_af + a(2*j-1) * a(2*j)
1224 ENDDO
1225 IF ( sum_mu_x_af <= 0.0) THEN
1226 IF (idebug > 0) THEN
1227 WRITE(iout,*) ' SUM_MU_X_AF = ', sum_mu_x_af
1228 ENDIF
1229 ivalid = 0
1230 GOTO 49
1231 ENDIF
1232 ELSEIF (icheck == 2) THEN
1233 DO j = 1, nmual/2
1234 IF ( ( a(2*j-1) * a(2*j) ) <= 0.0 ) THEN
1235 IF (idebug > 0) THEN
1236 WRITE(iout,*) ' A(2*J-1) * A(2*J) <= 0 '
1237 WRITE(iout,*) ' J, A(2*J-1), A(2*J) = ',
1238 $ j, a(2*j-1), a(2*j)
1239 ENDIF
1240 ivalid = 0
1241 GOTO 49
1242 ENDIF
1243 ENDDO
1244 ENDIF
1245
124649 CONTINUE
1247
1248 IF (idebug > 0) THEN
1249 WRITE(iout,*) 'IVALID = ', ivalid
1250 ENDIF
1251
1252 RETURN
1253
1254 END SUBROUTINE law69_check
1255
1256!> @purpose
1257!> Esstimate lwer/upper bounds of mu_i and alpha_i
1258!> which will be used to generate initial guess
1259!> @param[in] LAWID 1=Ogden, 2=Mooney-Rivlin
1260!> @param[in] ICHECK Checking level
1261!> @param[in] X(NPT) stretch values of given stress/strain(stretch) curve
1262!> @param[in] Y(NPT) engineering stress values of given stress/strain curve
1263!> @param[in] NPT number of points on the given stress/strain curve
1264
1265!> @param[out] MCOF_MIN(NMUAL) estimated lower bounds of A(1:NMUAL)
1266!> @param[out] MCOF_MAX(NMUAL) estimated upper bounds of A(1:NMUAL)
1267!||====================================================================
1268!|| law69_guess_bounds ../starter/source/materials/tools/nlsqf.F
1269!||--- called by ------------------------------------------------------
1270!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
1271!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
1272!||====================================================================
1273 SUBROUTINE law69_guess_bounds(LAWID, ICHECK, X, Y, NPT,
1274 $ NMUAL, MCOF_MIN, MCOF_MAX,ICOMP)
1275#include "implicit_f.inc"
1276#include "units_c.inc"
1277
1278 INTEGER LAWID, ICHECK, NPT,ICOMP
1279 my_real X(NPT), Y(NPT)
1280 INTEGER NMUAL
1281 my_real MCOF_MIN(NMUAL), MCOF_MAX(NMUAL)
1282
1283 INTEGER I
1284
1285 my_real AVE_SLOPE, mu_max, dx,DY
1286
1287 ave_slope = zero
1288 DO i = 1, npt
1289 dx = x(i) - one
1290c avolid dx to be too small
1291 IF (dx >= zero) THEN
1292 dx = max(dx, em6)
1293 ELSE
1294 dx = min(dx, -em6)
1295 ENDIF
1296 ave_slope = ave_slope + (y(i) - zero) / dx
1297!ok AVE_SLOPE = MAX(AVE_SLOPE , (Y(I) - ZERO) / dx )
1298 ENDDO
1299 ave_slope = ave_slope / (one * npt)
1300 IF (idebug > 0) THEN
1301 WRITE(iout, *) 'AVE_SLOPE = ', ave_slope
1302 ENDIF
1303 mu_max = ave_slope
1304 IF(icomp == 0 ) mu_max = max(ave_slope, twenty)
1305 IF (lawid == 1) THEN
1306 DO i = 1, nmual/2
1307 ! mu_i
1308 mcof_min(2*i-1) = -mu_max
1309 mcof_max(2*i-1) = mu_max
1310
1311 ! alpha_i
1312 mcof_min(2*i) = -10.0
1313 mcof_max(2*i) = +10.0
1314 ENDDO
1315 ELSEIF (lawid == 2) THEN
1316
1317 IF (icheck == 2) THEN
1318 mcof_min(1) = 0.0
1319 mcof_max(1) = mu_max
1320 ELSE
1321 mcof_min(1) = -mu_max
1322 mcof_max(1) = mu_max
1323 ENDIF
1324
1325 mcof_min(2) = +2.0
1326 mcof_max(2) = +2.0
1327
1328 IF (icheck == 2) THEN
1329 mcof_min(3) = -mu_max
1330 mcof_max(3) = 0.0
1331 ELSE
1332 mcof_min(3) = -mu_max
1333 mcof_max(3) = mu_max
1334 ENDIF
1335
1336 mcof_min(4) = -2.0
1337 mcof_max(4) = -2.0
1338
1339 ENDIF
1340
1341 IF (idebug > 0) THEN
1342 DO i = 1, nmual
1343 WRITE(iout,*) 'I, MCOF_MIN(I), MCOF_MAX(I) = ',
1344 $ i, mcof_min(i), mcof_max(i)
1345 ENDDO
1346 ENDIF
1347
1348 END SUBROUTINE law69_guess_bounds
1349
1350!> @purpose
1351!> return a random value within [0.0, 1.0) with uniform distribution
1352!> @param[in] IVAL integer for calculating next random number
1353!> @note
1354!> - random number generator of Park and Miller.
1355!> - Schrage's modification is used to avoid 'a*(m-1)' overflows for 32-bit.
1356!||====================================================================
1357!|| rand1 ../starter/source/materials/tools/nlsqf.F
1358!||--- called by ------------------------------------------------------
1359!|| law69_guess1 ../starter/source/materials/tools/nlsqf.f
1360!||--- calls -----------------------------------------------------
1361!|| my_exit ../starter/source/output/analyse/analyse.c
1362!||====================================================================
1363 REAL function rand1(ival)
1364#include "implicit_f.inc"
1365 INTEGER IVAL
1366 INTEGER IA, IM, IQ, IR
1367 REAL RIM
1368 PARAMETER (IA=69621, IM=2147483647,RIM=1./IM,IQ=30845,IR=23902)
1369
1370 INTEGER IDQ
1371
1372 IF (IVAL <= 0) THEN
1373 WRITE(*,*) 'ERROR, IVAL <= 0'
1374 WRITE(*,*) __FILE__,__LINE__
1375 CALL MY_EXIT(2)
1376 ELSEIF (IVAL > IM) THEN
1377 WRITE(*,*) 'ERROR, IVAL > IM'
1378 WRITE(*,*) __FILE__,__LINE__
1379 CALL MY_EXIT(2)
1380 ENDIF
1381
1382 IDQ = IVAL / IQ
1383 IVAL = IA * (IVAL - IDQ * IQ) - IR * IDQ
1384 if (IVAL < 0) IVAL = IVAL + IM
1385
1386 RAND1 = REAL(IVAL) * RIM
1387
1388 RETURN
1389 END FUNCTION RAND1
my_exit
void my_exit(int *i)
Definition analyse.c:1038
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
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
zero2small
subroutine zero2small(value, small)
Definition nlsqf.F:1079
ogden1
subroutine ogden1(x, a, dyda, ia, na)
@purpose calculate d(stress)/d(A) OGDEN (uniaxial tension)
Definition nlsqf.F:508
rand1
real function rand1(ival)
@purpose return a random value within [0.0, 1.0) with uniform distribution
Definition nlsqf.F:1364
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_nlsqf
subroutine law69_nlsqf(lawid, stretch, y, ma, nmual, npt, mual, icheck, nstart, errtol, id, titr, is_encrypted)
Definition nlsqf.F:565
inversion
subroutine inversion(a, n, np, b, m, mp, iret)
Definition nlsqf.F:304
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
mrqcof
subroutine mrqcof(x, y, sig, ndata, a, ia, ma, alpha, beta, nalp, errnow, ifuncs, icomp)
Definition nlsqf.F:188
covsrt
subroutine covsrt(covar, sizcovar, ma, ia, mfit)
Definition nlsqf.F:420
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
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