OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
nl_solv.F File Reference
#include "implicit_f.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"
#include "impl1_c.inc"
#include "impl2_c.inc"
#include "units_c.inc"
#include "task_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine nl_solv (nddl, iddl, ndof, ikc, d, dr, nnz, iadk, jdik, diag_k, lt_k, f, nddli, iadi, jdii, diag_i, lt_i, itok, iadm, jdim, diag_m, lt_m, r02, dd, ddr, itask0, it, itc, ru0, rold, idiv, inprint, icprec, istop, e02, de0, eimp, inloc, nddl0, ls, u02, gap, itab, fr_elem, iad_elem, w_ddl, a, ar, v, ms, x, ipari, intbuf_tab, num_imp, ns_imp, ne_imp, nsrem, nsl, icont, graphe, fac_k, ipiv_k, nk, nmonv, imonv, monvol, igrsurf, fr_mv, volmon, ibfv, skew, xframe, mumps_par, cddlp, ind_imp, nbintc, intlist, newfront, isendto, irecvfrom, irbe3, lrbe3, ndiv, icont0, isign, fext, dg, dgr, dg0, dgr0, rfext, ls1, nodft, nodlt, irbe2, lrbe2, idiv0, relres, anew_stif)
subroutine crit_ite (it, ur, rr, er, ndiv, tol)
subroutine line_s (r0, s0, r1, s, dtol, idiv, iint, prec)
subroutine get_max (n, d, vmax, i)
subroutine line_s1 (e1, s, ep, sp, en, sn, idiv, dtol, icont, icont0, iriks)
subroutine al_constraint1 (nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2, ier)
subroutine al_constraint2 (nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2)
subroutine b_correct (nddl, f, fext, dlamda)
subroutine b_correct_h (f_ddl, l_ddl, f, fext, dlamda)
subroutine al_constrainth1 (nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2, ft_n, lt_n, f_ddl, l_ddl, itask, ier)
subroutine al_constrainth2 (nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2, ft_n, lt_n, f_ddl, l_ddl, itask)
subroutine b_correct_hp (nddl, f, fext, dlamda)
subroutine al_constraint1_hp (nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2, ier)
subroutine al_constraint2_hp (nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2)

Function/Subroutine Documentation

◆ al_constraint1()

subroutine al_constraint1 ( integer nddl0,
integer nddl,
integer, dimension(*) iddl,
integer, dimension(*) ndof,
integer, dimension(*) ikc,
dd,
ddr,
dg,
dgr,
di,
dir,
integer, dimension(*) w_ddl,
l_a,
lamda,
sw2,
integer ier )

Definition at line 1116 of file nl_solv.F.

1120C-----------------------------------------------
1121C I m p l i c i t T y p e s
1122C-----------------------------------------------
1123#include "implicit_f.inc"
1124C-----------------------------------------------
1125C C o m m o n B l o c k s
1126C-----------------------------------------------
1127#include "com04_c.inc"
1128C-----------------------------------------------
1129C D u m m y A r g u m e n t s
1130C-----------------------------------------------
1131 INTEGER NDDL,NDDL0,IDDL(*) ,NDOF(*) ,IKC(*),W_DDL(*) ,IER
1132C REAL
1133 my_real
1134 . lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
1135 . l_a,sw2
1136C-----------------------------------------------
1137C L o c a l V a r i a b l e s
1138C-----------------------------------------------
1139 INTEGER I
1140 my_real
1141 . a,b,c,ui2,uiud,uiut,ut2,lam1,lam2,s,th1,th2,fac
1142C-----------------------------------------------
1143 ier=0
1144 CALL produt_u(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1145 . dg ,dgr ,ut2 ,w_ddl )
1146 a = ut2+sw2
1147C-------DD< DI+DD
1148 CALL integrator2(1,numnod,iddl,ndof,ikc,dd,ddr,di,dir)
1149 CALL produt_u2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1150 . dg ,dgr ,dd ,ddr ,uiut ,
1151 . w_ddl )
1152 b= two*(uiut+sw2*lamda)
1153 CALL produt_u(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1154 . dd ,ddr ,ui2 ,w_ddl )
1155 c = ui2-l_a*l_a+sw2*lamda
1156 s= b*b-four*a*c
1157 IF (s>=0) THEN
1158 s=sqrt(s)
1159 lam1=half*(-b+s)/a
1160 lam2=half*(-b-s)/a
1161 CALL produt_u2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1162 . dd ,ddr ,di ,dir ,uiud ,
1163 . w_ddl )
1164 CALL produt_u2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1165 . dg ,dgr ,di ,dir ,uiut ,
1166 . w_ddl )
1167 th1= uiud+lam1*uiut
1168 th2= uiud+lam2*uiut
1169 IF (th1>zero.AND.th2<zero) THEN
1170 lamda =lam1
1171 ELSEIF (th1<zero.AND.th2>zero) THEN
1172 lamda =lam2
1173 ELSEIF (th1>zero.AND.th2>zero) THEN
1174 lamda =min(lam1,lam2)
1175 ELSE
1176 ier=1
1177 ENDIF
1178 ELSE
1179 ier=1
1180 ENDIF
1181 fac =-one
1182 CALL frac_dd(1,numnod,iddl,ndof,ikc,dd,ddr,di,dir,fac)
1183C
1184 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
integrator2
subroutine integrator2(nodft, nodlt, iddl, ndof, ikc, d, dr, dd, ddr)
Definition integrator.F:207
frac_dd
subroutine frac_dd(nodft, nodlt, iddl, ndof, ikc, d, dr, dd, ddr, fac)
Definition integrator.F:295
min
#define min(a, b)
Definition macros.h:20
produt_u
subroutine produt_u(nddl0, nddl, iddl, ndof, ikc, dd, ddr, norm2, w_imp)
Definition produt_v.F:224
produt_u2
subroutine produt_u2(nddl0, nddl, iddl, ndof, ikc, d1, d1r, d2, d2r, norm2, w_imp)
Definition produt_v.F:262

◆ al_constraint1_hp()

subroutine al_constraint1_hp ( integer nddl0,
integer nddl,
integer, dimension(*) iddl,
integer, dimension(*) ndof,
integer, dimension(*) ikc,
dd,
ddr,
dg,
dgr,
di,
dir,
integer, dimension(*) w_ddl,
l_a,
lamda,
sw2,
integer ier )

Definition at line 1460 of file nl_solv.F.

1464C-----------------------------------------------
1465C I m p l i c i t T y p e s
1466C-----------------------------------------------
1467#include "implicit_f.inc"
1468C-----------------------------------------------
1469C D u m m y A r g u m e n t s
1470C-----------------------------------------------
1471 INTEGER NDDL,NDDL0,IDDL(*) ,NDOF(*) ,IKC(*),W_DDL(*) ,IER,
1472 . FT_N ,LT_N
1473C REAL
1474 my_real
1475 . lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
1476 . l_a,sw2
1477C-----------------------------------------------
1478C L o c a l V a r i a b l e s
1479C-----------------------------------------------
1480 INTEGER I
1481 my_real
1482 . a,b,c,ui2,uiud,uiut,ut2,lam1,lam2,s,th1,th2,fac
1483C-----------------------------------------------
1484 ier=0
1485 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1486 . dg ,dgr ,ut2 ,w_ddl )
1487C---------------------
1488 a = ut2+sw2
1489C-------DD< DI+DD
1490 CALL integrator2_hp(iddl,ndof,ikc,dd,ddr,di,dir)
1491C---------------------
1492 CALL produt_uhp2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1493 . dg ,dgr ,dd ,ddr ,uiut ,
1494 . w_ddl )
1495 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1496 . dd ,ddr ,ui2 ,w_ddl )
1497C---------------------
1498 b= two*(uiut+sw2*lamda)
1499 c = ui2-l_a*l_a+sw2*lamda
1500 s= b*b-four*a*c
1501 IF (s>=0) THEN
1502 s=sqrt(s)
1503 lam1=half*(-b+s)/a
1504 lam2=half*(-b-s)/a
1505 CALL produt_uhp2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1506 . dd ,ddr ,di ,dir ,uiud ,
1507 . w_ddl )
1508 CALL produt_uhp2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1509 . dg ,dgr ,di ,dir ,uiut ,
1510 . w_ddl )
1511C---------------------
1512 th1= uiud+lam1*uiut
1513 th2= uiud+lam2*uiut
1514 IF (th1>zero.AND.th2<zero) THEN
1515 lamda =lam1
1516 ELSEIF (th1<zero.AND.th2>zero) THEN
1517 lamda =lam2
1518 ELSEIF (th1>zero.AND.th2>zero) THEN
1519 lamda =min(lam1,lam2)
1520 ELSE
1521 ier=1
1522 ENDIF
1523 ELSE
1524 ier=1
1525 ENDIF
1526 fac =-one
1527 CALL frac_dd_hp(iddl,ndof,ikc,dd,ddr,di,dir,fac)
1528C
1529 RETURN
frac_dd_hp
subroutine frac_dd_hp(iddl, ndof, ikc, d, dr, dd, ddr, fac)
Definition integrator.F:618
integrator2_hp
subroutine integrator2_hp(iddl, ndof, ikc, d, dr, dd, ddr)
Definition integrator.F:513
produt_uhp
subroutine produt_uhp(nddl0, nddl, iddl, ndof, ikc, dd, ddr, norm2, w_imp)
Definition produt_v.F:3359
produt_uhp2
subroutine produt_uhp2(nddl0, nddl, iddl, ndof, ikc, d1, d1r, d2, d2r, norm2, w_imp)
Definition produt_v.F:3398

◆ al_constraint2()

subroutine al_constraint2 ( integer nddl0,
integer nddl,
integer, dimension(*) iddl,
integer, dimension(*) ndof,
integer, dimension(*) ikc,
dd,
ddr,
dg,
dgr,
di,
dir,
integer, dimension(*) w_ddl,
l_a,
lamda,
sw2 )

Definition at line 1192 of file nl_solv.F.

1195C-----------------------------------------------
1196C I m p l i c i t T y p e s
1197C-----------------------------------------------
1198#include "implicit_f.inc"
1199C-----------------------------------------------
1200C D u m m y A r g u m e n t s
1201C-----------------------------------------------
1202 INTEGER NDDL,NDDL0,IDDL(*) ,NDOF(*) ,IKC(*),W_DDL(*)
1203 my_real
1204 . lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
1205 . l_a,sw2
1206C-----------------------------------------------
1207C L o c a l V a r i a b l e s
1208C-----------------------------------------------
1209 INTEGER I
1210 my_real
1211 . gi,ui2,uiud,uiut,s
1212C-----------------------------------------------
1213 CALL produt_u(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1214 . di ,dir ,ui2 ,w_ddl )
1215 gi = sqrt(ui2)
1216 CALL produt_u2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1217 . dg ,dgr ,di ,dir ,uiut ,
1218 . w_ddl )
1219 CALL produt_u2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1220 . dd ,ddr ,di ,dir ,uiud ,
1221 . w_ddl )
1222 s = sw2*lamda
1223 lamda = ((l_a-gi)*gi-uiud)/(uiut+s)
1224C
1225 RETURN

◆ al_constraint2_hp()

subroutine al_constraint2_hp ( integer nddl0,
integer nddl,
integer, dimension(*) iddl,
integer, dimension(*) ndof,
integer, dimension(*) ikc,
dd,
ddr,
dg,
dgr,
di,
dir,
integer, dimension(*) w_ddl,
l_a,
lamda,
sw2 )

Definition at line 1539 of file nl_solv.F.

1542C-----------------------------------------------
1543C I m p l i c i t T y p e s
1544C-----------------------------------------------
1545#include "implicit_f.inc"
1546C-----------------------------------------------
1547C D u m m y A r g u m e n t s
1548C-----------------------------------------------
1549 INTEGER NDDL,NDDL0,IDDL(*) ,NDOF(*) ,IKC(*),W_DDL(*)
1550C REAL
1551 my_real
1552 . lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
1553 . l_a,sw2
1554C-----------------------------------------------
1555C L o c a l V a r i a b l e s
1556C-----------------------------------------------
1557 INTEGER I
1558 my_real
1559 . gi,ui2,uiud,uiut,s
1560C-----------------------------------------------
1561 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1562 . di ,dir ,ui2 ,w_ddl )
1563 CALL produt_uhp2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1564 . dg ,dgr ,di ,dir ,uiut ,
1565 . w_ddl )
1566 CALL produt_uhp2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1567 . dd ,ddr ,di ,dir ,uiud ,
1568 . w_ddl )
1569C---------------------
1570 gi = sqrt(ui2)
1571 s = sw2*lamda
1572 lamda = ((l_a-gi)*gi-uiud)/(uiut+s)
1573C
1574 RETURN

◆ al_constrainth1()

subroutine al_constrainth1 ( integer nddl0,
integer nddl,
integer, dimension(*) iddl,
integer, dimension(*) ndof,
integer, dimension(*) ikc,
dd,
ddr,
dg,
dgr,
di,
dir,
integer, dimension(*) w_ddl,
l_a,
lamda,
sw2,
integer ft_n,
integer lt_n,
integer f_ddl,
integer l_ddl,
integer itask,
integer ier )

Definition at line 1287 of file nl_solv.F.

1292C-----------------------------------------------
1293C I m p l i c i t T y p e s
1294C-----------------------------------------------
1295#include "implicit_f.inc"
1296C-----------------------------------------------
1297C D u m m y A r g u m e n t s
1298C-----------------------------------------------
1299 INTEGER NDDL,NDDL0,IDDL(*) ,NDOF(*) ,IKC(*),W_DDL(*) ,IER,
1300 . FT_N ,LT_N ,F_DDL ,L_DDL ,ITASK
1301C REAL
1302 my_real
1303 . lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
1304 . l_a,sw2
1305C-----------------------------------------------
1306C L o c a l V a r i a b l e s
1307C-----------------------------------------------
1308 INTEGER I
1309 my_real
1310 . a,b,c,ui2,uiud,uiut,ut2,lam1,lam2,s,th1,th2,fac
1311C-----------------------------------------------
1312 ier=0
1313 CALL produt_uh(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1314 . dg ,dgr ,ut2 ,w_ddl ,f_ddl ,
1315 . l_ddl ,itask )
1316C----------------------
1317 CALL my_barrier
1318C---------------------
1319 a = ut2+sw2
1320C-------DD< DI+DD
1321 CALL integrator2(ft_n,lt_n,iddl,ndof,ikc,dd,ddr,di,dir)
1322C----------------------
1323 CALL my_barrier
1324C---------------------
1325 CALL produt_uh2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1326 . dg ,dgr ,dd ,ddr ,uiut ,
1327 . w_ddl ,f_ddl ,l_ddl ,itask )
1328 CALL produt_uh(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1329 . dd ,ddr ,ui2 ,w_ddl ,f_ddl ,
1330 . l_ddl ,itask )
1331C----------------------
1332 CALL my_barrier
1333C---------------------
1334 b= two*(uiut+sw2*lamda)
1335 c = ui2-l_a*l_a+sw2*lamda
1336 s= b*b-four*a*c
1337 IF (s>=0) THEN
1338 s=sqrt(s)
1339 lam1=half*(-b+s)/a
1340 lam2=half*(-b-s)/a
1341 CALL produt_uh2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1342 . dd ,ddr ,di ,dir ,uiud ,
1343 . w_ddl ,f_ddl ,l_ddl ,itask )
1344 CALL produt_uh2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1345 . dg ,dgr ,di ,dir ,uiut ,
1346 . w_ddl ,f_ddl ,l_ddl ,itask )
1347C----------------------
1348 CALL my_barrier
1349C---------------------
1350 th1= uiud+lam1*uiut
1351 th2= uiud+lam2*uiut
1352 IF (th1>zero.AND.th2<zero) THEN
1353 lamda =lam1
1354 ELSEIF (th1<zero.AND.th2>zero) THEN
1355 lamda =lam2
1356 ELSEIF (th1>zero.AND.th2>zero) THEN
1357 lamda =min(lam1,lam2)
1358 ELSE
1359 ier=1
1360 ENDIF
1361 ELSE
1362 ier=1
1363 ENDIF
1364 fac =-one
1365 CALL frac_dd(ft_n,lt_n,iddl,ndof,ikc,dd,ddr,di,dir,fac)
1366C
1367 RETURN
produt_uh2
subroutine produt_uh2(nddl0, nddl, iddl, ndof, ikc, d1, d1r, d2, d2r, norm2, w_imp, f_ddl, l_ddl, itask)
Definition produt_v.F:1741
produt_uh
subroutine produt_uh(nddl0, nddl, iddl, ndof, ikc, dd, ddr, norm2, w_imp, f_ddl, l_ddl, itask)
Definition produt_v.F:1687
my_barrier
subroutine my_barrier
Definition machine.F:31

◆ al_constrainth2()

subroutine al_constrainth2 ( integer nddl0,
integer nddl,
integer, dimension(*) iddl,
integer, dimension(*) ndof,
integer, dimension(*) ikc,
dd,
ddr,
dg,
dgr,
di,
dir,
integer, dimension(*) w_ddl,
l_a,
lamda,
sw2,
integer ft_n,
integer lt_n,
integer f_ddl,
integer l_ddl,
integer itask )

Definition at line 1376 of file nl_solv.F.

1380C-----------------------------------------------
1381C I m p l i c i t T y p e s
1382C-----------------------------------------------
1383#include "implicit_f.inc"
1384C-----------------------------------------------
1385C D u m m y A r g u m e n t s
1386C-----------------------------------------------
1387 INTEGER NDDL,NDDL0,IDDL(*) ,NDOF(*) ,IKC(*),W_DDL(*),
1388 . FT_N ,LT_N ,F_DDL ,L_DDL ,ITASK
1389C REAL
1390 my_real
1391 . lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
1392 . l_a,sw2
1393C-----------------------------------------------
1394C L o c a l V a r i a b l e s
1395C-----------------------------------------------
1396 INTEGER I
1397 my_real
1398 . gi,ui2,uiud,uiut,s
1399C-----------------------------------------------
1400 CALL produt_uh(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1401 . di ,dir ,ui2 ,w_ddl ,f_ddl ,
1402 . l_ddl ,itask )
1403 CALL produt_uh2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1404 . dg ,dgr ,di ,dir ,uiut ,
1405 . w_ddl ,f_ddl ,l_ddl ,itask )
1406 CALL produt_uh2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
1407 . dd ,ddr ,di ,dir ,uiud ,
1408 . w_ddl ,f_ddl ,l_ddl ,itask )
1409C----------------------
1410 CALL my_barrier
1411C---------------------
1412 gi = sqrt(ui2)
1413 s = sw2*lamda
1414 lamda = ((l_a-gi)*gi-uiud)/(uiut+s)
1415C
1416 RETURN

◆ b_correct()

subroutine b_correct ( integer nddl,
f,
fext,
dlamda )

Definition at line 1230 of file nl_solv.F.

1231C-----------------------------------------------
1232C I m p l i c i t T y p e s
1233C-----------------------------------------------
1234#include "implicit_f.inc"
1235C-----------------------------------------------
1236C D u m m y A r g u m e n t s
1237C-----------------------------------------------
1238 INTEGER NDDL
1239 my_real
1240 . dlamda,f(*),fext(*)
1241C-----------------------------------------------
1242C L o c a l V a r i a b l e s
1243C-----------------------------------------------
1244 INTEGER I
1245C-----------------------------------------------
1246 DO i =1,nddl
1247 f(i)=f(i)+dlamda*fext(i)
1248 ENDDO
1249C
1250 RETURN

◆ b_correct_h()

subroutine b_correct_h ( integer f_ddl,
integer l_ddl,
f,
fext,
dlamda )

Definition at line 1255 of file nl_solv.F.

1256C-----------------------------------------------
1257C I m p l i c i t T y p e s
1258C-----------------------------------------------
1259#include "implicit_f.inc"
1260C-----------------------------------------------
1261C D u m m y A r g u m e n t s
1262C-----------------------------------------------
1263 INTEGER F_DDL ,L_DDL
1264C REAL
1265 my_real
1266 . dlamda,f(*),fext(*)
1267C-----------------------------------------------
1268C L o c a l V a r i a b l e s
1269C-----------------------------------------------
1270 INTEGER I
1271C-----------------------------------------------
1272 DO i =f_ddl ,l_ddl
1273 f(i)=f(i)+dlamda*fext(i)
1274 ENDDO
1275C
1276 RETURN

◆ b_correct_hp()

subroutine b_correct_hp ( integer nddl,
f,
fext,
dlamda )

Definition at line 1423 of file nl_solv.F.

1424C-----------------------------------------------
1425C I m p l i c i t T y p e s
1426C-----------------------------------------------
1427#include "implicit_f.inc"
1428C-----------------------------------------------
1429C D u m m y A r g u m e n t s
1430C-----------------------------------------------
1431 INTEGER NDDL
1432C REAL
1433 my_real
1434 . dlamda,f(*),fext(*)
1435C-----------------------------------------------
1436C L o c a l V a r i a b l e s
1437C-----------------------------------------------
1438 INTEGER F_DDL ,L_DDL ,ITSK
1439 INTEGER I
1440C-----------------------------------------------
1441!$OMP PARALLEL PRIVATE(ITSK,F_DDL ,L_DDL,I)
1442 CALL imp_smpini(itsk ,f_ddl ,l_ddl ,nddl )
1443 DO i =f_ddl ,l_ddl
1444 f(i)=f(i)+dlamda*fext(i)
1445 ENDDO
1446!$OMP END PARALLEL
1447C
1448 RETURN
imp_smpini
subroutine imp_smpini(itsk, n1ftsk, n1ltsk, n1)
Definition imp_solv.F:6910

◆ crit_ite()

subroutine crit_ite ( integer it,
ur,
rr,
er,
integer ndiv,
tol )

Definition at line 809 of file nl_solv.F.

810C-----------------------------------------------
811C I m p l i c i t T y p e s
812C-----------------------------------------------
813#include "implicit_f.inc"
814C-----------------------------------------------
815C C o m m o n B l o c k s
816C-----------------------------------------------
817#include "impl1_c.inc"
818#include "impl2_c.inc"
819C-----------------------------------------------
820C D u m m y A r g u m e n t s
821C-----------------------------------------------
822 INTEGER IT,NDIV
823C REAL
824 my_real
825 . ur,rr,er,tol
826C-----------------------------------------------
827C L o c a l V a r i a b l e s
828C-----------------------------------------------
829 my_real
830 . err,tolr,told
831C--------cas particulier---------------------------------------
832C TOLD=EP04
833C IF (ILINE_S==1) TOLD=EP03
834C IF (ISMDISP==1) TOLD=EP10
835 told = tol_div
836 tolr = one+tol
837C--------1:E,2:r,3:E+U,4:r+U,5,r+E,6,r+E+U-----
838C--new: 1:E,2:r,3:U,12:E+r,13,E+U,23,r+U,123,r+E+U
839 IF (nitol==1) THEN
840 err=abs(er)/tol
841 IF (it==1.OR.rr>=one) err=err+tolr
842 ELSEIF (nitol==2) THEN
843 err=rr/tol
844 ELSEIF (nitol==3) THEN
845 err=ur/tol
846C ERR=MAX(SQRT(ABS(ER)),UR)
847 IF (rr>=one) err=err+tolr
848 ELSEIF (nitol==12) THEN
849 err=max(rr/n_tolf,abs(er)/n_tole)
850c ERR=MAX(RR,UR)
851 ELSEIF (nitol==13) THEN
852 err=max(ur/n_tolu,abs(er)/n_tole)
853c ERR=MAX(SQRT(ABS(ER)),RR)
854 ELSEIF (nitol==23) THEN
855 err=max(ur/n_tolu,rr/n_tolf)
856c ERR=MAX(SQRT(ABS(ER)),RR,UR)
857 ELSEIF (nitol==123) THEN
858 err=max(ur/n_tolu,rr/n_tolf,abs(er)/n_tole)
859 ENDIF
860C
861 IF (err<=one) THEN
862C IF (ERR<=TOL) THEN
863 imconv=1
864!sb ELSEIF (UR<MIN(EM05,TOL).AND.RR<EM01) THEN
865 ELSEIF ((n_lim /= 1.OR.isprb /= 0).AND.ur<min(em05,tol).AND.rr<em01) THEN
866 imconv=1
867 ELSEIF (imconv<0) THEN
868C--------line-search-----
869 ELSE
870 imconv=0
871C------divergence-----------
872!sb
873 IF(n_lim /= 1 .OR. isprb /= 0) THEN
874 IF (rr>tolr) THEN
875 IF (it==1) THEN
876 imconv=-3
877 IF (ndiv<ndiver) imconv=0
878 ELSE
879 imconv=-2
880 IF ((ndiv<ndiver.OR.ilast==1).AND.rr<told) imconv=0
881 ENDIF
882 ENDIF
883 ENDIF
884
885 IF (idtc>0.AND.it>nl_dtn) imconv=-2
886 ENDIF
887C-----------------------------------------------
888 RETURN
max
#define max(a, b)
Definition macros.h:21

◆ get_max()

subroutine get_max ( integer n,
d,
vmax,
integer i )

Definition at line 962 of file nl_solv.F.

963C-----------------------------------------------
964C I m p l i c i t T y p e s
965C-----------------------------------------------
966#include "implicit_f.inc"
967C-----------------------------------------------
968C C o m m o n B l o c k s
969C-----------------------------------------------
970#include "com01_c.inc"
971C-----------------------------------------------
972C D u m m y A r g u m e n t s
973C-----------------------------------------------
974 INTEGER N,I
975 my_real
976 . d(3,*),vmax
977C-----------------------------------------------
978C L o c a l V a r i a b l e s
979C-----------------------------------------------
980 INTEGER K
981 my_real
982 . t
983C-----------------------------------------------
984 vmax=zero
985 i=1
986 DO k=1,n
987 t=max(abs(d(1,k)),abs(d(2,k)),abs(d(3,k)))
988 IF (t>vmax) THEN
989 vmax = t
990 i=k
991 ENDIF
992 ENDDO
993 IF (nspmd>1) CALL spmd_max_s(vmax)
994C
995 RETURN
spmd_max_s
subroutine spmd_max_s(s)
Definition imp_spmd.F:1195

◆ line_s()

subroutine line_s ( r0,
s0,
r1,
s,
dtol,
integer idiv,
integer iint,
prec )

Definition at line 895 of file nl_solv.F.

896C-----------------------------------------------
897C I m p l i c i t T y p e s
898C-----------------------------------------------
899#include "implicit_f.inc"
900C-----------------------------------------------
901C C o m m o n B l o c k s
902C-----------------------------------------------
903#include "impl1_c.inc"
904C-----------------------------------------------
905C D u m m y A r g u m e n t s
906C-----------------------------------------------
907 INTEGER IDIV,IINT
908C REAL
909 my_real
910 . r0,r1,s0,s,dtol,prec
911C-----------------------------------------------
912C L o c a l V a r i a b l e s
913C-----------------------------------------------
914 INTEGER I
915 my_real
916 . s1,dr,r
917C-----------------------------------------------
918C R1=ABS(E1/E0)
919C IF (R1<=ONE) THEN
920 dr=abs(r0-r1)
921 r=min(r1/prec,dr)
922c R=MIN(R1,DR*PREC)
923 IF ((r<dtol.OR.idiv<-nls_lim)
924 . .OR.(r1>one.AND.iint>0.AND.irefi<2)) THEN
925 IF (imconv==-1)imconv=0
926c IDIV=0
927 idiv=idiv-1
928 RETURN
929 ELSEIF (r1>r0.AND.imconv==-1) THEN
930 idiv=-nls_lim-2
931 s=s0
932 RETURN
933 ENDIF
934 imconv=-1
935 s1=s
936 s=abs(s0*r1-s1*r0)/dr
937 IF (r1>one) THEN
938C----------diverge---------------------------
939 IF (s>one) s=half*s1/r1
940 s=max(em01,s)
941C IDIV=IDIV+1
942 ELSE
943C----- extrapolation --------
944 IF (s>max(s1,s0).OR.s<min(s1,s0)) THEN
945 s=min(three_half*s1,s)
946 s=min(five,s)
947 ENDIF
948 ENDIF
949 s0=s1
950 r0=r1
951 idiv=idiv-1
952C
953 RETURN

◆ line_s1()

subroutine line_s1 ( e1,
s,
ep,
sp,
en,
sn,
integer idiv,
dtol,
integer icont,
integer icont0,
integer iriks )

Definition at line 1002 of file nl_solv.F.

1003C-----------------------------------------------
1004C I m p l i c i t T y p e s
1005C-----------------------------------------------
1006#include "implicit_f.inc"
1007C-----------------------------------------------
1008C C o m m o n B l o c k s
1009C-----------------------------------------------
1010#include "impl1_c.inc"
1011C-----------------------------------------------
1012C D u m m y A r g u m e n t s
1013C-----------------------------------------------
1014 INTEGER IDIV,ICONT,ICONT0,IRIKS
1015C REAL
1016 my_real
1017 . e1,s,ep,sp,en,sn,dtol
1018C-----------------------------------------------
1019C L o c a l V a r i a b l e s
1020C-----------------------------------------------
1021 INTEGER I,IM
1022 my_real
1023 . sold,stmp,s2,tol,eold
1024C-------min e=r.du-----------------E1-> already relative /E0-----------------------
1025C------find sl so that el<0
1026 imconv =-1
1027 eold=zero
1028 im = max(icont,icont0)
1029 IF (e1>zero) THEN
1030 IF (idiv==0) THEN
1031 idiv = idiv + 1
1032 sp = one
1033 ep = e1
1034 s = s + one
1035 IF (iriks>0.OR.im > 0) THEN
1036 s = one
1037 imconv =0
1038 END IF
1039 RETURN
1040 ELSEIF (idiv>0) THEN
1041 sold = s
1042 IF (abs(e1-ep)<dtol.OR.e1>ep.OR.idiv>nls_lim) THEN
1043 IF (e1>ep) s = s - one
1044 imconv =0
1045 ELSE
1046 s = s + one
1047 ENDIF
1048 idiv = idiv + 1
1049 sp = sold
1050 ep = e1
1051 RETURN
1052 ENDIF
1053 sp = s
1054 ep = e1
1055 ELSE
1056 IF (idiv==0) THEN
1057 sp = zero
1058 ep = one
1059 s = half
1060 idiv = idiv - 1
1061 sn = one
1062 en = e1
1063 RETURN
1064 ELSEIF (sp==zero) THEN
1065 IF (abs(e1-en)<dtol.OR.idiv<-nls_lim) THEN
1066 imconv =0
1067 ELSEIF (iriks>0.OR.ismdisp>0) THEN
1068 imconv =0
1069 ELSE
1070c SN = S
1071c EN = E1
1072 s = s*half
1073 ENDIF
1074 idiv = idiv - 1
1075 RETURN
1076 ELSEIF (idiv>0) THEN
1077 idiv = -idiv
1078 ENDIF
1079 eold=en
1080 sn = s
1081 en = e1
1082 ENDIF
1083 idiv = idiv - 1
1084c----------return to
1085 sold = s
1086 im = max(icont,icont0)
1087 IF (im>0) THEN
1088c S2 = (EP-EN)/SN+SN
1089c STMP = HALF*(S2-SQRT(S2*S2-FOUR*EP))
1090 stmp = half*(-sp+sn)
1091 ELSE
1092 stmp = ep*(sn-sp)/(ep-en)
1093 ENDIF
1094 s = sp + stmp
1095 stmp =(sold-s)/sold
1096 stmp = max(abs(stmp),abs(e1-eold))
1097 IF (stmp<dtol.OR.idiv<-nls_lim) THEN
1098 imconv =0
1099 ELSEIF (iriks>0.AND.
1100 . (s>two.OR.s<em01)) THEN
1101 imconv =0
1102 ELSE
1103 imconv =-1
1104 ENDIF
1105C
1106 RETURN

◆ nl_solv()

subroutine nl_solv ( integer nddl,
integer, dimension(*) iddl,
integer, dimension(*) ndof,
integer, dimension(*) ikc,
d,
dr,
integer nnz,
integer, dimension(*) iadk,
integer, dimension(*) jdik,
diag_k,
lt_k,
f,
integer nddli,
integer, dimension(*) iadi,
integer, dimension(*) jdii,
diag_i,
lt_i,
integer, dimension(*) itok,
integer, dimension(*) iadm,
integer, dimension(*) jdim,
diag_m,
lt_m,
r02,
dd,
ddr,
integer itask0,
integer it,
integer itc,
ru0,
rold,
integer idiv,
integer inprint,
integer icprec,
integer istop,
e02,
de0,
eimp,
integer, dimension(*) inloc,
integer nddl0,
ls,
u02,
gap,
integer, dimension(*) itab,
integer, dimension(*) fr_elem,
integer, dimension(*) iad_elem,
integer, dimension(*) w_ddl,
a,
ar,
v,
ms,
x,
integer, dimension(*) ipari,
type(intbuf_struct_), dimension(*) intbuf_tab,
integer, dimension(*) num_imp,
integer, dimension(*) ns_imp,
integer, dimension(*) ne_imp,
integer nsrem,
integer nsl,
integer icont,
type(prgraph), dimension(*) graphe,
fac_k,
integer, dimension(*) ipiv_k,
integer nk,
integer nmonv,
integer, dimension(*) imonv,
integer, dimension(*) monvol,
type (surf_), dimension(nsurf) igrsurf,
integer, dimension(*) fr_mv,
volmon,
integer, dimension(*) ibfv,
skew,
xframe,
integer mumps_par,
integer, dimension(*) cddlp,
integer, dimension(*) ind_imp,
integer nbintc,
integer, dimension(*) intlist,
integer, dimension(*) newfront,
integer, dimension(*) isendto,
integer, dimension(*) irecvfrom,
integer, dimension(*) irbe3,
integer, dimension(*) lrbe3,
integer ndiv,
integer icont0,
integer isign,
fext,
dg,
dgr,
dg0,
dgr0,
rfext,
ls1,
integer nodft,
integer nodlt,
integer, dimension(*) irbe2,
integer, dimension(*) lrbe2,
integer idiv0,
relres,
character*1 anew_stif )

Definition at line 53 of file nl_solv.F.

74C-----------------------------------------------
75C M o d u l e s
76C-----------------------------------------------
77 USE dsgraph_mod
78 USE intbufdef_mod
79 USE groupdef_mod
80C-----------------------------------------------
81C I m p l i c i t T y p e s
82C-----------------------------------------------
83#include "implicit_f.inc"
84C-----------------------------------------------
85C C o m m o n B l o c k s
86C-----------------------------------------------
87#if defined(MUMPS5)
88#include "dmumps_struc.h"
89#endif
90#include "com01_c.inc"
91#include "com04_c.inc"
92#include "com08_c.inc"
93#include "impl1_c.inc"
94#include "impl2_c.inc"
95#include "units_c.inc"
96#include "task_c.inc"
97C-----------------------------------------------
98C D u m m y A r g u m e n t s
99C-----------------------------------------------
100C----------resol [K]{X}={F}------
101 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),IADM(*),JDIM(*),ITASK0,
102 . NDDLI ,IADI(*),JDII(*),ITOK(*),ICPREC,INLOC(*),
103 . NDOF(*),IDDL(*),IKC(*),IDIV ,INPRINT,ISTOP,
104 . IT,ITC,NDDL0,ITAB(*),FR_ELEM(*),IAD_ELEM(*),W_DDL(*)
105 INTEGER NE_IMP(*),NSREM ,NSL,ICONT,ISIGN,IMP1,
106 . IPARI(*) ,NUM_IMP(*),NS_IMP(*), IPIV_K(*), NK
107 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*),
108 . IBFV(*),IND_IMP(*),IRBE3(*) ,LRBE3(*),NDIV ,ICONT0,
109 . IRBE2(*),LRBE2(*)
110 INTEGER NEWFRONT(*),NBINTC,INTLIST(*),ISENDTO(*),IRECVFROM(*),
111 . NODFT ,NODLT,IDIV0
112C REAL
113 my_real
114 . diag_k(*),lt_k(*),f(*),r02,diag_m(*),lt_m(*),
115 . diag_i(*),lt_i(*),d(3,*),dr(3,*),dd(3,*),ddr(3,*),
116 . ru0,rold,e02 ,eimp,lstol,de0,ls(*),u02,gap,rbid
117 my_real
118 . a(3,*),ar(3,*),v(3,*),x(3,*),ms(*),fac_k(*),
119 . volmon(*),skew(*),xframe(*),fext(*),dg(3,*),dgr(3,*),
120 . dg0(3,*),dgr0(3,*),ls1(*) ,rfext, relres(*)
121C-----LS(1):R0,LS(2):S0,LS(3):S1 ; LS1(1): EN, LS1(2):SN;LS1(3):DE
122 TYPE(PRGRAPH) :: GRAPHE(*)
123 INTEGER CDDLP(*)
124 TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
125 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
126C
127#ifdef MUMPS5
128 TYPE(DMUMPS_STRUC) MUMPS_PAR
129#else
130 ! Fake declaration as DMUMPS_STRUC is shipped with MUMPS
131 INTEGER MUMPS_PAR
132#endif
133C
134 character*1 anew_stif
135#ifdef MUMPS5
136C-----------------------------------------------
137C L o c a l V a r i a b l e s
138C-----------------------------------------------
139C----------resol [K]{X}={F}------
140 INTEGER I ,J ,ND,IP,ITMP,ICOV,ICALU,IER,ISETK0,IRIKS,
141 . F_DDL ,L_DDL
142C REAL
143 my_real
144 . r2,u2,du2,tol,rr,ru,temp,fr,re,de,r1,nl_tol,r01,rr1,
145 . gapn,lamda,fa,uc2,um1,um2,pmax
146C-----------------------------------------------
147C------------------------------------------
148c insolv=1 => (elast) Newton modified
149c insolv=2,3 => BFGS
150c insolv=4,5 => elastoplastic
151C------------------------------------------
152 fr = 0 ! was not always initialized before use
153C------------------------------------------
154C------------------------------------------
155C [K]:stiffness matrix [M]:preconditioning matrix
156c F_DDL=1+ITASK*NDDL/NTHREAD
157c L_ddl = (itask+1)*nddl/nthread
158C
159 tol =l_tol
160 nl_tol=n_tol
161C
162 lstol=ls_tol
163 ip=iabs(inprint)
164 de=0
165 nd=3*numnod
166 ru=ru0
167 isetk0=isetk+icont
168 IF (nitol==3.OR.nitol==13.OR.nitol==23
169 . .OR.nitol==123) THEN
170 icalu=1
171 ELSE
172 icalu=0
173 ENDIF
174 iriks=0
175 IF (ncycle>1.AND.idtc==3) THEN
176 IF (ilast==0) iriks=1
177 ENDIF
178C----------------------
179c CALL MY_BARRIER
180C---------------------
181 isetk=0
182 IF (icont>0) gapn=gap*zep9
183 IF (imconv>=0) icont0=icont
184C--------idtc=3 :modified Riks arc_length-----
185 IF (it==0) THEN
186 r01=sqrt(r02)
187C--------ROLD=R2(IT=0)---------------------------------------
188Ctmp IF (ISPRB==1) R01=MIN(R01,SQRT(ROLD)/R01)
189 re=one
190C--------cas particulier---------------------------------------
191 IF (r01*nl_tol<=em12.AND.nl_tol>em10.AND.idiv/=10) THEN
192 IF (irefi/=5.OR.icont==0) THEN
193 IF(inprint/=0)THEN
194 WRITE(iout,*)
195c WRITE(IOUT,1008)IT,RE,R01,RE
196 WRITE(iout,1108)it,anew_stif,re,r01,re
197 WRITE(iout,*)
198 IF(inprint<0)THEN
199 WRITE(istdo,*)
200c WRITE(ISTDO,1008)IT,RE,R01,RE
201 WRITE(istdo,1108)it,anew_stif,re,r01,re
202 WRITE(istdo,*)
203 ENDIF
204 ENDIF
205 relres(1) = re
206 relres(2) = r01
207 relres(3) = re
208 idiv = -2
209 isetk=1
210 RETURN
211 ELSE
212 f(1)=max(em20,f(1))
213 r01=max(r01,em20)
214 ENDIF
215 ENDIF
216 IF (insolv==4) icov = imconv
217 imconv=0
218 rr=rold
219 ndiv=0
220C---------
221 IF (insolv==2.OR.insolv==3) CALL bfgs_0
222C---------
223 IF (ncycle==1) isign = 1
224 IF (iriks>0) THEN
225 CALL lin_solv(nddl ,iddl ,ndof ,ikc ,dg ,
226 1 dgr ,tol ,nnz ,iadk ,jdik ,
227 2 diag_k,lt_k ,nddli ,iadi ,jdii ,
228 3 diag_i,lt_i ,itok ,iadm ,jdim ,
229 4 diag_m,lt_m ,fext ,de ,inloc ,
230 5 fr_elem,iad_elem,w_ddl,itask0,icprec,
231 6 istop ,a ,ar ,v ,
232 7 ms ,x ,ipari ,intbuf_tab ,
233 8 num_imp,ns_imp,ne_imp,nsrem ,nsl ,
234 9 itc ,graphe, itab, fac_k, ipiv_k,
235 a nk ,nmonv,imonv ,monvol,igrsurf,
236 b fr_mv ,volmon,ibfv ,skew ,
237 c xframe ,mumps_par,cddlp,ind_imp,rbid,
238 e irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
239
240C---------
241 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
242 . dg ,dgr ,u2 ,w_ddl )
243C-------distinquer converge or diverge---
244citask0 IF (ITASK == 0) THEN
245 IF (idiv/=-2) THEN
246 CALL cp_real(nd,dd,dg0)
247 IF (iroddl/=0) CALL cp_real(nd,ddr,dgr0)
248 ELSE
249 dla_riks = isign*dt2
250 ENDIF
251citask0 END IF !(ITASK == 0) THEN
252 CALL produt_uhp2(nddl0 ,nddl ,iddl ,ndof ,ikc ,
253 . dg ,dgr ,dg0 ,dgr0 ,uc2 ,
254 . w_ddl )
255C----------------------
256c CALL MY_BARRIER
257C---------------------
258citask0 IF (ITASK == 0) THEN
259 CALL get_max(numnod,dg0,temp,i)
260 um1=temp
261 CALL get_max(numnod,dgr0,temp,i)
262 um1=max(um1,temp,em20)
263 CALL get_max(numnod,dg,temp,i)
264 um2=temp
265 CALL get_max(numnod,dgr,temp,i)
266 um2=max(um2,temp,em20)
267 uc2=uc2/um1/um2
268C LAMDA = ALEN/SQRT(U2)
269 IF ((uc2+dla_riks/rfext)<zero) THEN
270 isign = -1
271 ELSE
272 isign = 1
273 ENDIF
274c print *,'UC2,DLA_RIKS,ISIGN=',UC2,DLA_RIKS,ISIGN
275 IF (impdeb>0) THEN
276 IF (ncycle>=ndeb0.AND.ncycle<=ndeb1)
277 . WRITE(iout,1025)uc2,dla_riks,isign
278 ENDIF
279 IF (isign<0) THEN
280 dla_riks = 2*isign*dt2
281 ELSE
282 dla_riks = zero
283 ENDIF
284C---------
285citask0 END IF !(ITASK == 0) THEN
286
287C----------------------
288c CALL MY_BARRIER
289C---------------------
290 lamda = isign*dt2
291 IF (idiv/=-2.AND.icont==0)THEN
292 CALL produt_vmhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
293 . dd ,ddr ,f ,de ,w_ddl )
294 CALL integrator2_hp(iddl,ndof,ikc,d,dr,dd,ddr)
295 ELSE
296 CALL frac_dd_hp(iddl,ndof,ikc,d,dr,dg,dgr,lamda)
297 CALL produt_vmhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
298 . d ,dr ,f ,de ,w_ddl )
299 ENDIF
300 ELSE
301
302 IF ((n_lim /= 1.OR.isprb /= 0).AND.
303 . (ncycle>1.AND.idiv/=-2.AND.icont==0.AND.ikt==0
304 . .OR.(ncycle==1.AND.idiv==10)))THEN
305 IF (idtc==3.AND.ilast==1) THEN
306 lamda = dla_riks
307 dla_riks = zero
308 CALL frac_d_hp(iddl,ndof,ikc,dd,ddr,lamda)
309 ENDIF
310 CALL produt_vmhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
311 . dd ,ddr ,f ,de ,w_ddl )
312 ELSE
313 IF (insolv==4.AND.icov>=0) imconv=-1
314
315 CALL lin_solv(nddl ,iddl ,ndof ,ikc ,dd ,
316 1 ddr ,tol ,nnz ,iadk ,jdik ,
317 2 diag_k,lt_k ,nddli ,iadi ,jdii ,
318 3 diag_i,lt_i ,itok ,iadm ,jdim ,
319 4 diag_m,lt_m ,f ,de ,inloc ,
320 5 fr_elem,iad_elem,w_ddl,itask0,icprec,
321 6 istop ,a ,ar ,v ,
322 7 ms ,x ,ipari ,intbuf_tab ,
323 8 num_imp,ns_imp,ne_imp,nsrem ,nsl ,
324 9 itc ,graphe, itab, fac_k, ipiv_k,
325 a nk ,nmonv,imonv ,monvol,igrsurf,
326 b fr_mv ,volmon,ibfv ,skew ,
327 c xframe ,mumps_par,cddlp,ind_imp,rbid,
328 e irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
329 IF (de<zero.AND.irefi>1.AND.irefi<5) imconv=-2
330 IF (icont>0.AND.idsgap>0) THEN
331 CALL get_max(numnod,dd,pmax,i)
332C----------------------
333c CALL MY_BARRIER
334C---------------------
335 IF (pmax>gapn) THEN
336 temp=gapn/pmax
337 CALL frac_d_hp(iddl,ndof,ikc,dd,ddr,temp)
338 de = de * temp
339 ENDIF
340 ENDIF
341 ENDIF
342C
343 CALL integrator2_hp(iddl,ndof,ikc,d,dr,dd,ddr)
344 ENDIF !IF (IRIKS>0)
345C----------------------
346c CALL MY_BARRIER
347C---------------------
348citask0 IF (ITASK == 0) THEN
349 IF (imconv>=0) THEN
350 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
351 . d ,dr ,u2 ,w_ddl )
352 eimp=eimp+de
353 IF (de>ep10.AND.iline==0.AND.idyna==0
354 . .AND.ncycle==1) THEN
355 IF(ispmd==0)THEN
356 WRITE(iout,1030)eimp
357 WRITE(istdo,1030)eimp
358 ENDIF
359 CALL imp_stop(-1)
360 ENDIF
361C--------for idtc>=2-----
362 u02=u2
363 IF(inprint/=0)THEN
364 ru=sqrt(u2)
365c WRITE(IOUT,*)
366c WRITE(IOUT,1002)IT,RU,R01,EIMP
367 WRITE(iout,1102)
368 WRITE(iout,1003)it,anew_stif,ru,r01,eimp
369 IF(inprint<0)THEN
370c WRITE(ISTDO,*)
371c WRITE(ISTDO,1002)IT,RU,R01,EIMP
372 WRITE(istdo,1102)
373 WRITE(istdo,1003)it,anew_stif,ru,r01,eimp
374 ENDIF
375 ENDIF
376
377 relres(1) = ru
378 relres(2) = r01
379 relres(3) = eimp
380 ELSE
381 WRITE(iout,*)
382 WRITE(iout,1013) de
383 IF(inprint<0)THEN
384 WRITE(istdo,*)
385 WRITE(istdo,1013)de
386 ENDIF
387C
388 ENDIF !IF (IMCONV>=0)
389citask0 END IF !(ITASK == 0) THEN
390C---------IT>0---------
391 ELSE
392C--------R=Fext-Fint correction due to Fext---------
393 IF (iriks>0.AND.dla_riks/=zero) THEN
394 CALL vaxpy_hp(nddl, f ,fext,dla_riks)
395 ENDIF
396 CALL produt_hp(nddl,f,f,w_ddl,r2)
397C----------------------
398c CALL MY_BARRIER
399C---------------------
400 IF (r2>=zero.AND.r2<ep30) THEN
401 rr=sqrt(r2/r02)
402 IF (it==1.AND.isprb==1.AND.idiv/=-2) THEN
403 rr1=sqrt(r2/rold)
404C
405 IF (icont>0.AND.rr<one) rr1=rr
406 ELSE
407 rr1=rr
408 ENDIF
409 re=de0/eimp
410 IF (it==1) THEN
411 ru=one
412 ELSEIF (icalu>0) THEN
413 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
414 . d ,dr ,u2 ,w_ddl )
415 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
416 . dd ,ddr ,du2 ,w_ddl )
417C----------------------
418c CALL MY_BARRIER
419C---------------------
420 ru=ls(3)*sqrt(du2/u2)
421 IF (insolv==5) ru=sqrt(du2/u2)
422 ENDIF
423 CALL crit_ite(it,ru,rr1,re,ndiv,nl_tol)
424 ELSE
425 imconv =-2
426 rr1 =ep30
427 ENDIF
428C----------------------
429c CALL MY_BARRIER
430C---------------------
431 IF (imconv==-3) THEN
432 IF (ncycle==1.OR.idiv==-2) imconv=-2
433 r02=r2
434 ENDIF
435C----------------------
436c CALL MY_BARRIER
437C---------------------
438 IF (imconv==0) THEN
439 idiv=0
440C-------start line-search------------
441 IF (iline_s/=1) THEN
442 ls(1)=one
443 ls(2)=zero
444 ENDIF
445 ls(3)=one
446 IF (ndiver>0) r02=r2
447 ENDIF
448 IF(imconv==1) THEN
449 istop=0
450 idiv=0
451C------For time correction-with Riks----
452 IF (it>1.AND.icalu==0) THEN
453 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
454 . d ,dr ,u2 ,w_ddl )
455 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
456 . dd ,ddr ,du2 ,w_ddl )
457C----------------------
458c CALL MY_BARRIER
459C---------------------
460 ru=ls(3)*sqrt(du2/u2)
461 u02=u2
462 ELSEIF (idtc>=2) THEN
463 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
464 . d ,dr ,u02 ,w_ddl )
465C----------------------
466c CALL MY_BARRIER
467C---------------------
468 ENDIF
469C
470 IF(inprint/=0)THEN
471c WRITE(IOUT,*)
472c WRITE(IOUT,1008)IT,RU,RR,RE
473 WRITE(iout,1108)it,anew_stif,ru,rr,re
474 WRITE(iout,*)
475 IF(inprint<0)THEN
476c WRITE(ISTDO,*)
477c WRITE(ISTDO,1008)IT,RU,RR,RE
478 WRITE(istdo,1108)it,anew_stif,ru,rr,re
479 WRITE(istdo,*)
480 ENDIF
481 ENDIF
482 relres(1) = ru
483 relres(2) = rr
484 relres(3) = re
485 ELSEIF(imconv<=-2) THEN
486C--------diverge----
487C write(iout,*)'diverge',IMCONV,IT ,IDIV
488 IF(inprint/=0)THEN
489 WRITE(iout,*)
490 IF (rr1>one) THEN
491 WRITE(iout,1011)rr1
492 ELSEIF(it>nl_dtn) THEN
493 WRITE(iout,1010)
494 ELSEIF(idtc==3) THEN
495 WRITE(iout,1020)
496 ENDIF
497 WRITE(iout,*)
498 IF(inprint<0)THEN
499 WRITE(istdo,*)
500 IF (rr1>one) THEN
501 WRITE(istdo,1011)rr1
502 ELSEIF(it>nl_dtn) THEN
503 WRITE(istdo,1010)
504 ELSEIF(idtc==3) THEN
505 WRITE(istdo,1020)
506 ENDIF
507 WRITE(istdo,*)
508 ENDIF
509 IF(imconv==-2) THEN
510 WRITE(iout,1012)
511 IF(inprint<0)WRITE(istdo,1012)
512 ENDIF
513 ENDIF
514 relres(1) = ru
515 relres(2) = rr1
516 relres(3) = re
517 ELSE
518C--------CAS IMCONV=0,-1---add line-search for constant dt-
519c IF (IDTC>0) THEN
520 IF (iline_s>0.AND.isign>=0.AND.irwall==0) THEN
521 fr=ls(3)
522 IF (nitol/=2.AND.nitol/=4.OR.iline_s==1) THEN
523C----------DD,DDR are modified in some cases
524 IF (ls1(3)/=zero ) THEN
525 de = ls1(3)
526 ELSE
527 CALL produt_vmhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
528 . dd ,ddr ,f ,de ,w_ddl )
529 END IF
530C----------------------
531c CALL MY_BARRIER
532C---------------------
533 r1=de/de0
534 IF (iline_s==3) THEN
535 r1=abs(de/de0)
536 IF (rr>one) r1=one
537 IF (rr>one.AND.irefi>=2) r1=ls(1)
538 temp=ep02
539 ENDIF
540 ELSE
541 r1=rr/rold
542 temp=ep03
543 ENDIF
544 icov = imconv
545C------------
546citask0 IF (ITASK == 0) THEN
547C------------
548 IF (iline_s==3) THEN
549 CALL line_s(ls(1),ls(2),r1,fr,lstol,idiv,icont,temp)
550 ELSEIF (iline_s==1) THEN
551 CALL line_s1(r1,fr,ls(1),ls(2),ls1(1),ls1(2),idiv,lstol,
552 . icont,icont0,iriks)
553 ELSEIF (iline_s==2) THEN
554 CALL line_s(ls(1),ls(2),r1,fr,lstol,idiv,icont,temp)
555 ENDIF
556 IF (impdeb>0.AND.ispmd==0) THEN
557 IF (ncycle>=ndeb0.AND.ncycle<=ndeb1) then
558 WRITE(iout,*)'R1,FR,IDIV=',r1,fr,idiv
559C WRITE(IOUT,*)'DE,DE0=',DE,DE0
560 end if
561 ENDIF
562C------------
563citask0 END IF !(ITASK == 0) THEN
564C----------------------
565c CALL MY_BARRIER
566C---------------------
567C
568 IF ((insolv==2.OR.insolv==3).AND.fr>one.AND.ikt==0)
569 . fr = min(fr,two)
570 ENDIF
571C IF (IABS(INSOLV)==5) THEN
572C----------suprimed----------
573C ELSE
574 IF (imconv==0) THEN
575C----------schema normale---------
576 IF (it>1) THEN
577C------------
578 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
579 . d ,dr ,u2 ,w_ddl )
580 CALL produt_uhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
581 . dd ,ddr ,du2 ,w_ddl )
582C----------------------
583c CALL MY_BARRIER
584C---------------------
585 ru=ls(3)*sqrt(du2/u2)
586 ENDIF
587C------------
588citask0 IF (ITASK == 0) THEN
589C------------
590 IF (insolv==2.OR.insolv==3) CALL bfgs_ls(ls(3))
591 IF (rr>=rold) THEN
592 ndiv = ndiv +1
593 ELSE
594 ndiv = 0
595 ENDIF
596 IF (inprint/=0)THEN
597 IF(mod(it,ip)==0)THEN
598c WRITE(IOUT,*)
599 WRITE(iout,1003)it,anew_stif,ru,rr,re
600 IF(inprint<0)THEN
601c WRITE(ISTDO,*)
602 WRITE(istdo,1003)it,anew_stif,ru,rr,re
603 ENDIF
604 ENDIF
605 ENDIF
606 relres(1) = ru
607 relres(2) = rr
608 relres(3) = re
609 IF(itc>=n_lim.AND.ismdisp==0) THEN
610c IF(INPRINT/=0)THEN
611c WRITE(IOUT,*)
612c WRITE(IOUT,1006)ITC
613c IF(INPRINT<0)THEN
614c WRITE(ISTDO,*)
615c WRITE(ISTDO,1006)ITC
616c ENDIF
617c ENDIF
618 itc=0
619 ENDIF
620C------------
621citask0 END IF !(ITASK == 0) THEN
622C----------------------
623c CALL MY_BARRIER
624C---------------------
625 IF (iriks>0) THEN
626 IF (isetk0==1) THEN
627 CALL lin_solv(nddl ,iddl ,ndof ,ikc ,dg ,
628 1 dgr ,tol ,nnz ,iadk ,jdik ,
629 2 diag_k,lt_k ,nddli ,iadi ,jdii ,
630 3 diag_i,lt_i ,itok ,iadm ,jdim ,
631 4 diag_m,lt_m ,fext ,de ,inloc ,
632 5 fr_elem,iad_elem,w_ddl,itask0,icprec,
633 6 istop ,a ,ar ,v ,
634 7 ms ,x ,ipari ,intbuf_tab ,
635 8 num_imp,ns_imp,ne_imp,nsrem ,nsl ,
636 9 itc ,graphe, itab, fac_k, ipiv_k,
637 a nk ,nmonv,imonv ,monvol,igrsurf,
638 b fr_mv ,volmon,ibfv ,skew ,
639 c xframe ,mumps_par,cddlp,ind_imp,rbid,
640 e irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
641 icprec = 0
642 idsc = 0
643 END IF
644 CALL produt_vmhp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
645 . dg ,dgr ,f ,um1 ,w_ddl )
646 END IF
647C
648 CALL lin_solv(nddl ,iddl ,ndof ,ikc ,dd ,
649 1 ddr ,tol ,nnz ,iadk ,jdik ,
650 2 diag_k,lt_k ,nddli ,iadi ,jdii ,
651 3 diag_i,lt_i ,itok ,iadm ,jdim ,
652 4 diag_m,lt_m ,f ,de ,inloc ,
653 5 fr_elem,iad_elem,w_ddl,itask0,icprec,
654 6 istop ,a ,ar ,v ,
655 7 ms ,x ,ipari ,intbuf_tab ,
656 8 num_imp,ns_imp,ne_imp,nsrem ,nsl ,
657 9 itc ,graphe, itab, fac_k, ipiv_k,
658 a nk ,nmonv,imonv ,monvol,igrsurf,
659 b fr_mv ,volmon,ibfv ,skew ,
660 c xframe ,mumps_par,cddlp,ind_imp,rbid,
661 d irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
662 IF (icont>0.AND.idsgap>0) THEN
663 CALL get_max(numnod,dd,pmax,i)
664C----------------------
665c CALL MY_BARRIER
666C---------------------
667 IF (pmax>gapn) THEN
668 temp=gapn/pmax
669 CALL frac_d_hp(iddl,ndof,ikc,dd,ddr,temp)
670 de = de * temp
671 ENDIF
672 ENDIF
673C----------------------
674c CALL MY_BARRIER
675C---------------------
676C---------Riks type arc length------
677 IF (iriks>0) THEN
678 lamda = dla_riks
679 IF (ial_m==1) THEN
680 CALL al_constraint1_hp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
681 . dd ,ddr ,dg ,dgr ,d ,
682 . dr ,w_ddl ,alen ,lamda ,scal_riks,
683 . ier )
684C----------------------
685c CALL MY_BARRIER
686C---------------------
687 IF (ier==1) THEN
688 imconv=-2
689 IF(inprint/=0)THEN
690 WRITE(iout,*)
691 WRITE(iout,1020)
692 WRITE(iout,*)
693 IF(inprint<0)THEN
694 WRITE(istdo,*)
695 WRITE(istdo,1020)
696 WRITE(istdo,*)
697 ENDIF
698 ENDIF
699 END IF !(IER==1) THEN
700 ELSEIF (ial_m==2) THEN
701 CALL al_constraint2_hp(nddl0 ,nddl ,iddl ,ndof ,ikc ,
702 . dd ,ddr ,dg ,dgr ,d ,
703 . dr ,w_ddl ,alen ,lamda ,scal_riks)
704C----------------------
705c CALL MY_BARRIER
706C---------------------
707 END IF
708 IF (imconv>=0) THEN
709 CALL frac_dd_hp(iddl,ndof,ikc,d ,dr ,dg,dgr,lamda)
710citask0 IF (ITASK == 0) THEN
711 dla_riks = dla_riks + lamda
712 eimp=eimp+dla_riks*um1
713citask0 END IF !(ITASK == 0) THEN
714C----------------------
715c CALL MY_BARRIER
716C---------------------
717 ENDIF
718 END IF !IF (IRIKS>0) THEN
719C
720 CALL integrator2_hp(iddl,ndof,ikc,d,dr,dd,ddr)
721 eimp=eimp+de
722c if (de<0) print *,'****negative de=',de
723C
724 ELSEIF(imconv==-1) THEN
725C----------line-search---------
726 temp=-ls(3)+fr
727 CALL frac_dd_hp(iddl,ndof,ikc,d,dr,dd,ddr,temp)
728C----------------------
729c CALL MY_BARRIER
730C---------------------
731 ls(3) = fr
732 END IF !IF (IMCONV==0) THEN
733C
734 END IF !IF(IMCONV==1) THEN
735C-----fin if IT=0
736 ENDIF
737 IF (imconv>=0) THEN
738 ru0=ru
739 de0=de
740 IF (abs(de0)<em20) de0=em20
741 rold=max(em20,rr)
742 ENDIF
743C-----REFORMING [K]----
744 IF (ismdisp>0) THEN
745 IF (itc>=n_lim.AND.imconv>=0) THEN
746 itc = 0
747 isetk=1
748 END IF
749 IF (itc==0.AND.imconv>=0) THEN
750 IF (ncycle==1.OR.n_lim==1) isetk=1
751 IF (isolv==5.OR.isolv==6) isetk=1
752 END IF
753 IF (irwall >0 .AND.imconv == 1) isetk=1
754 IF ((itc==0.OR.imconv==1).AND.ikt>0) isetk=1
755 ELSE
756C----condition IMCONV==1 reforming [K] for nothing (no resolution)-
757 IF (itc==0.OR.imconv==1) isetk=1
758 END IF
759 IF ((idtc==3.OR.ikt>0).AND.imconv<=-2) isetk=1
760C-----if still diverge, try w/o [K] reforming
761 IF (imconv<=-2.AND.rr1>ep20.AND.idiv0>=0) isetk=1
762 IF ((istop==1.OR.istop==2).AND.(isolv==5.OR.isolv==6))THEN
763 istop = 0
764 isetk = 1
765 imconv=-2
766 ENDIF
767 IF (imconv<=-2.AND.isprb>0.AND.rr1>ep04) isetk=1
768! sb
769 IF(imconv<= -2.AND.n_lim == 1 .AND. isprb == 0) isetk = 1
770
771C----------------------
772c CALL MY_BARRIER
773C---------------------
774Ctmp IF (ITC==0.OR.IMCONV==1.OR.IMCONV<=-2) ISETK=1
775c 1002 FORMAT(5X,'NL_ITERATION=',I5,1X,' INITIAL RESIDUAL NORM=',3E11.4)
776 1102 FORMAT(3x,78('-')/
777 . 12x,'Stif. Mat.'/
778 . 6x,'Iter',2x, ' reformed ',5x,'|du|/|u|',3x,'|r|/|r0|',3x,'|dE|/|E|',1x,'Conv.stat.'/
779 . 3x,78('-'))
780c 1003 FORMAT(5X,'NL_ITERATION=',I5,1X,' RELATIVE RESIDUAL NORM=',3E11.4)
781 1003 FORMAT(5x,i5,6x,a1,5x,3(1x,1pe10.3))
782 1005 FORMAT(3x,'TOLERANCE FOR LINEAR ITERATIVE SOLVER :',2x,e11.4)
783 1006 FORMAT(3x,'--STIFFNESS MATRIX IS RESET AFTER ',i4,
784 . ' ITERATIONS--')
785 1007 FORMAT(3x,'--ITERATION DIVERGE with R=',e11.4,2x,
786 . 'STIFFNESS MATRIX WILL BE RESET')
787c 1008 FORMAT(3X,'* CONVERGED WITH ',I4,1X,'ITERATIONS,',
788c . ' |du|/|u|,|r|/|r0|,|dE|/|E|=',3E11.4)
789 1108 FORMAT(5x,i5,6x,a1,5x,3(1x,1pe10.3),5x, 'C')
790 1009 FORMAT(3x,'--ITERATION DIVERGE with RELATIVE R=',e11.4/,
791 . 3x,'RESET ITERATION WITH STABILIZATION BY SCALING= ',e11.4)
792 1010 FORMAT(3x,'--ITERATION DIVERGE with MAX_ITER REACHED--')
793 1011 FORMAT(3x,'ITERATION DIVERGE with RELATIVE R=',e11.4)
794 1012 FORMAT(3x,'--RESET ITERATION WITH NEW TIMESTEP--')
795 1013 FORMAT(3x,'ITERATION DIVERGE with NEGATIVE ENERGY DE=',e11.4)
796 1020 FORMAT(3x,'--ITERATION DIVERGE with RIKS Method--')
797 1025 FORMAT(3x,'UC2,DLA_RIKS,ISIGN=',2e11.4,i4)
798 1030 FORMAT(3x,'CHECK CONSTRAINT CONDITIONS, TOO LARGE ENERGY VALUE=',
799 . 2x,e11.4)
800 RETURN
801#endif
bfgs_0
subroutine bfgs_0
Definition imp_bfgs.F:81
bfgs_ls
subroutine bfgs_ls(ls)
Definition imp_bfgs.F:109
imp_stop
subroutine imp_stop(istop)
Definition imp_solv.F:1997
frac_d_hp
subroutine frac_d_hp(iddl, ndof, ikc, d, dr, fac)
Definition integrator.F:565
lin_solv
subroutine lin_solv(nddl, iddl, ndof, ikc, d, dr, tol, nnz, iadk, jdik, diag_k, lt_k, nddli, iadi, jdii, diag_i, lt_i, itok, iadm, jdim, diag_m, lt_m, f, f_u, inloc, fr_elem, iad_elem, w_ddl, itask, icprec, istop, a, ar, ve, ms, xe, ipari, intbuf_tab, num_imp, ns_imp, ne_imp, nsrem, nsl, it, graphe, itab, fac_k, ipiv_k, nk, nmonv, imonv, monvol, igrsurf, fr_mv, volmon, ibfv, skew, xframe, mumps_par, cddlp, ind_imp, xi_c, irbe3, lrbe3, irbe2, lrbe2)
Definition lin_solv.F:74
dsgraph_mod
Definition dsgraph_mod.F:31
groupdef_mod
Definition groupdef_mod.F:668
al_constraint2_hp
subroutine al_constraint2_hp(nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2)
Definition nl_solv.F:1542
crit_ite
subroutine crit_ite(it, ur, rr, er, ndiv, tol)
Definition nl_solv.F:810
al_constraint1_hp
subroutine al_constraint1_hp(nddl0, nddl, iddl, ndof, ikc, dd, ddr, dg, dgr, di, dir, w_ddl, l_a, lamda, sw2, ier)
Definition nl_solv.F:1464
get_max
subroutine get_max(n, d, vmax, i)
Definition nl_solv.F:963
line_s1
subroutine line_s1(e1, s, ep, sp, en, sn, idiv, dtol, icont, icont0, iriks)
Definition nl_solv.F:1003
line_s
subroutine line_s(r0, s0, r1, s, dtol, idiv, iint, prec)
Definition nl_solv.F:896
cp_real
subroutine cp_real(n, x, xc)
Definition produt_v.F:871
produt_hp
subroutine produt_hp(nddl, x, y, w, r)
Definition produt_v.F:3252
vaxpy_hp
subroutine vaxpy_hp(n, v, y, s)
Definition produt_v.F:3580
produt_vmhp
subroutine produt_vmhp(nddl0, nddl, iddl, ndof, ikc, dd, ddr, y, r, w_imp)
Definition produt_v.F:3321