00001 #include <math.h>
00002 #include <string.h>
00003 #include "pcpcsolver.h"
00004 #include "pglobalc.h"
00005 #include "globalc.h"
00006 #include "global.h"
00007 #include "globalt.h"
00008 #include "globmat.h"
00009 #include "globmatt.h"
00010 #include "globmatc.h"
00011 #include "elemswitcht.h"
00012 #include "gmatrix.h"
00013 #include "gtopology.h"
00014 #include "mechprint.h"
00015 #include "transprint.h"
00016 #include "backupsol.h"
00017 #include "backupsolt.h"
00018 #include "pcpsolver.h"
00019 #include "pglobal.h"
00020 #include "seqfilesm.h"
00021 #include "pglobalt.h"
00022 #include "seqfilest.h"
00023 #include "genfile.h"
00024 #include "mpi.h"
00025
00026
00027
00028
00029
00030
00031
00032
00033 void par_solve_gpcouplprob ()
00034 {
00035 long lcid;
00036
00037
00038
00039 lcid=0;
00040
00041 switch (Cp->tnlinsol){
00042 case newtonc:{
00043 par_newton_raphson_gparcoupl (lcid);
00044 break;
00045 }
00046 default:{
00047 print_err("unknown solver of nonlinear equation system is required",__FILE__,__LINE__,__func__);
00048 }
00049 }
00050 }
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 void par_newton_raphson_gparcoupl (long lcid)
00063 {
00064 switch (Tp->tprob){
00065 case growing_np_problem:{
00066
00067
00068 par_newton_raphson_gparcoupl_lin (lcid);
00069 break;
00070 }
00071 case growing_np_problem_nonlin:{
00072
00073
00074 par_newton_raphson_gparcoupl_nonlin (lcid);
00075 break;
00076 }
00077 default:{
00078 print_err("unknown TRFEL problem type is required",__FILE__,__LINE__,__func__);
00079 }
00080 }
00081 }
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 void par_newton_raphson_gparcoupl_lin (long lcid)
00094 {
00095 int mefel_convergence;
00096 long i,ii,j,k,nm,nt,ini,nsts,mncd,tncd,mnce,tnce,nii,mnae,nifn;
00097 long *ifn;
00098 double zero,dt,dtmin,dtdef,end_time,alpha,s,prev_timem,prev_timet,btime;
00099 double *dr,*fb,*r,*d,*f,*fp,*fi,*fl,*p,*lhst,*tdlhst,*rhst,*lhsb,*lhstb,*tdlhstb;
00100 double norf, norfa, err, lsm_res;
00101 matrix lsm_a(3,3);
00102 vector lsm_r(3), lsm_l(3);
00103
00104
00105 ini = Mp->niilnr;
00106
00107 err = Mp->errnr;
00108
00109 nm=Ndofm;
00110
00111 nt=Ndoft;
00112
00113
00114
00115
00116
00117
00118 lhst=Lsrst->give_lhs (lcid);
00119
00120 tdlhst=Lsrst->give_tdlhs (lcid);
00121
00122 rhst=Lsrst->give_rhs (lcid);
00123
00124
00125 d = new double [nt];
00126 nullv (d,nt);
00127
00128 p = new double [nt];
00129 nullv (p,nt);
00130
00131 lhstb = new double [nt];
00132 nullv (lhstb,nt);
00133
00134 tdlhstb = new double [nt];
00135 nullv (tdlhstb,nt);
00136
00137 ifn = new long[Mt->nn];
00138 memset(ifn, 0, sizeof(*ifn)*Mt->nn);
00139
00140
00141
00142
00143
00144
00145 r = Lsrs->give_lhs (lcid);
00146
00147 f = Lsrs->give_rhs (lcid);
00148
00149
00150 dr = new double [nm];
00151 nullv (dr,nm);
00152
00153 fp = new double [nm];
00154 nullv (fp,nm);
00155
00156 fl = new double [nm];
00157 nullv (fl,nm);
00158
00159 fi = new double [nm];
00160 nullv (fi,nm);
00161
00162 fb = new double [nm];
00163 nullv (fb,nm);
00164
00165 lhsb = new double [nm];
00166 nullv (lhsb,nm);
00167
00168
00169 alpha=Tp->alpha;
00170 zero=Tp->zero;
00171
00172
00173
00174 Cp->time=Cp->timecon.starttime ();
00175
00176 dt=Cp->timecon.initialtimeincr ();
00177
00178 end_time = Cp->timecon.endtime ();
00179
00180
00181 Mp->time = Cp->time;
00182 Mp->timecon.take_values (Cp->timecon);
00183
00184 dtmin=Mp->timecon.dtmin;
00185
00186 prev_timem = Mp->time;
00187
00188 Tp->time = Cp->time;
00189 Tp->timecont.take_values (Cp->timecon);
00190
00191 prev_timet = Tp->time;
00192
00193
00194 approximation_inittemper();
00195
00196
00197 approximation();
00198
00199 actual_previous_nodval();
00200
00201
00202
00203
00204 i=0;
00205 ii=0;
00206
00207 nsts=0;
00208
00209
00210 if (Mp->hdbcont.restore_stat()){
00211 if (Mespr==1)
00212 fprintf (stdout,"\n Reading of backup file for MEFEL\n");
00213 solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm);
00214 Mt->rhs_save (fp);
00215
00216 approximation_inittemper ();
00217
00218
00219
00220 Mm->initmaterialmodels(lcid);
00221 print_init(-1, "at",Pcp->fnim,Pcp->feim);
00222
00223
00224 }
00225 else
00226 {
00227
00228 approximation_inittemper ();
00229
00230
00231
00232 Mm->initmaterialmodels(lcid);
00233 print_init(-1, "wt",Pcp->fnim,Pcp->feim);
00234 print_step(lcid, i, Mp->time, f);
00235 print_flush();
00236 }
00237 if (Tp->hdbcont.restore_stat()){
00238 if (Mesprt==1)
00239 fprintf (stdout,"\n Reading of backup file for TRFEL\n");
00240 solvert_restore (lhst,tdlhst,f,i,Cp->time,dt,Cp->timecon,Ndoft);
00241
00242 Tm->initmaterialmodels ();
00243 compute_req_valt (0);
00244 print_initt(-1, "at",Pcp->fnit,Pcp->feit);
00245
00246
00247 }
00248 else{
00249
00250 Tm->initmaterialmodels ();
00251 compute_req_valt (0);
00252 print_initt(-1, "wt",Pcp->fnit,Pcp->feit);
00253 print_stept(0,i,Tp->time,NULL);
00254 print_flusht();
00255 }
00256
00257
00258 mefel_convergence = 1;
00259
00260
00261 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00262
00263 Mt->save_nodval(lcid);
00264
00265
00266
00267
00268
00269 do{
00270 if(mefel_convergence == 0){
00271 Cp->timecon.oldtime ();
00272 Tp->timecont.oldtime ();
00273 if (Mesprc==1) fprintf (stdout,"\n\n --------------------------------------------------------------------------------------------");
00274 if (Mesprc==1) fprintf (stdout,"\n REDUCED TIME STEP, METR TIME STEP = %ld : METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
00275 if (Mesprc==1) fprintf (stdout,"\n --------------------------------------------------------------------------------------------");
00276 i--;
00277 }
00278
00279
00280
00281 Cp->time = Cp->timecon.newtime (dt);
00282 Tp->time = Tp->timecont.newtime (dt);
00283 Tp->time = Cp->time;
00284 Tp->timecont.take_values (Cp->timecon);
00285 Tp->timecont = Cp->timecon;
00286 i++;
00287
00288 if(mefel_convergence == 1){
00289 if (Mesprc==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------------------");
00290 if (Mesprc==1) fprintf (stdout,"\n NEW METR TIME STEP = %ld: METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
00291 if (Mesprc==1) fprintf (stdout,"\n -------------------------------------------------------------------------------------------");
00292 }
00293 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
00294 if (Mesprt==1) fprintf (stdout,"\n TRFEL TIME STEP = %ld, TRFEL TIME = %e, trfel time increment = %e",i,Tp->time,dt);
00295 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
00296
00297
00298
00299
00300 tnce = Gtt->search_newelem (Tp->time,prev_timet);
00301
00302
00303
00304
00305 if (tnce>0){
00306 Tt->initial_nodval ();
00307 }
00308
00309
00310
00311
00312 tncd = Gtt->search_newdofs (Tp->time,prev_timet);
00313
00314
00315 Gtt->update_elem (Tp->time);
00316
00317
00318 Gtt->update_nodes ();
00319
00320
00321 Gtt->update_dofs (Tp->time);
00322
00323
00324
00325 nt=Ndoft;
00326
00327
00328 if (tncd>0){
00329 if (Kmat != NULL){
00330 delete Kmat;
00331 Kmat=NULL;
00332 }
00333
00334 if (Cmat != NULL){
00335 delete Cmat;
00336 Cmat=NULL;
00337 }
00338 }
00339
00340 if (tnce!=0 || tncd!=0){
00341
00342 Tt->lhs_restore (lhst,Lsrst->lhsi,tdlhst);
00343 }
00344
00345 Tm->updateipval ();
00346
00347
00348 conductivity_matrix (lcid);
00349
00350
00351 capacity_matrix (lcid);
00352
00353
00354 for (j=0;j<nt;j++){
00355 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
00356 }
00357
00358
00359 trfel_right_hand_side (0,rhst,nt);
00360
00361
00362 Kmat->gmxv (p,d);
00363
00364
00365
00366 Kmat->scalgm (dt*alpha);
00367 Kmat->addgm (1.0,*Cmat);
00368
00369 for (j=0;j<nt;j++){
00370 rhst[j] = rhst[j] - d[j];
00371 d[j]=tdlhst[j];
00372 }
00373
00374
00375 Psolt->par_linear_solver (Gtt,Kmat,tdlhst,rhst,Outt,Mesprt);
00376
00377
00378 for (j=0;j<nt;j++){
00379 lhstb[j]=lhst[j];
00380 tdlhstb[j]=tdlhst[j];
00381 }
00382
00383
00384 for (j=0;j<nt;j++){
00385 s=(1.0-alpha)*d[j]+alpha*tdlhst[j];
00386 lhst[j]+=dt*s;
00387 }
00388
00389
00390 solution_correction ();
00391
00392 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00393
00394
00395 approximation ();
00396
00397 actual_previous_nodval ();
00398
00399 if (Cp->fcsolv == fullnewtonc){
00400
00401
00402 if (Kmat != NULL){
00403 delete Kmat;
00404 Kmat=NULL;
00405 }
00406
00407 if (Cmat != NULL){
00408 delete Cmat;
00409 Cmat=NULL;
00410 }
00411 }
00412
00413
00414
00415
00416 mefel_convergence = 1;
00417
00418 Mp->time = Mp->timecon.newtime (dt);
00419
00420 if (Mp->time <= Tp->time){
00421
00422 ii++;
00423
00424
00425 Mp->strainstate=0;
00426
00427
00428 Mp->stressstate=0;
00429
00430
00431 Mp->otherstate=0;
00432
00433
00434 for (j=0;j<nm;j++){
00435 lhsb[j]=r[j];
00436 }
00437
00438 if (Mespr==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------");
00439 if (Mespr==1) fprintf (stdout,"\n MEFEL TIME STEP = %ld, MEFEL TIME = %e, mefel time increment = %e",ii,Mp->time,dt);
00440 if (Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------------------------");
00441
00442
00443
00444 mnce = Gtm->search_changed_elem(Mp->time, mnae);
00445
00446
00447 mncd = Gtm->search_changed_dofs(Mp->time,prev_timem);
00448
00449 nifn = 0;
00450 ifn = NULL;
00451 if (mnce!=0){
00452
00453
00454 nifn = Gtm->search_iface_nodes(ifn);
00455 }
00456
00457
00458 if ((mnce!=0) || (mncd!=0))
00459 {
00460
00461
00462 Mt->save_node_inidispl(lcid, Mp->time, prev_timem);
00463
00464
00465
00466
00467
00468 if (nifn)
00469 {
00470
00471 Gtm->switch_new_elem();
00472
00473
00474 Gtm->update_nodes();
00475
00476
00477
00478 Gtm->update_active_dofs(Mp->time);
00479
00480
00481 Ndofm = Gtm->codenum_generation(Out);
00482
00483
00484
00485 Mt->restore_nodval(r, Ndofm);
00486
00487
00488 stress_initdispl(lcid);
00489
00490
00491
00492 Gtm->update_active_dofs(Mp->time);
00493
00494 Gtm->clear_intf_dofs(ifn);
00495
00496
00497 Ndofm = Gtm->codenum_generation(Out);
00498
00499
00500 Mp->strcomp = 0;
00501 internal_forces(lcid,fb);
00502 cmulv(-1.0, fb, Ndofm);
00503 Mp->strcomp = 1;
00504
00505
00506
00507 delete Smat;
00508 Smat=NULL;
00509 stiffness_matrix (lcid);
00510 nullv(r, Ndofm);
00511
00512
00513
00514
00515 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00516
00517
00518 addv(r, dr, nm);
00519
00520
00521
00522 Mt->save_elem_inidispl(lcid, ifn);
00523
00524
00525
00526 Gtm->remove_nodes(ifn);
00527
00528
00529
00530 Mt->save_nodval(lcid);
00531 }
00532
00533
00534 Gtm->update_elem(Mp->time);
00535
00536
00537 Gtm->update_nodes();
00538
00539
00540 Gtm->update_dofs (Mp->time);
00541
00542
00543 nm = Ndofm = Gtm->codenum_generation(Out);
00544
00545
00546 delete Smat;
00547 Smat=NULL;
00548
00549
00550 Mt->restore_nodval(r,nm);
00551
00552
00553
00554
00555
00556
00557 Gtm->update_elem(prev_timem);
00558
00559 Gtm->update_nodes();
00560
00561 btime = Mp->time;
00562
00563 Mp->time = prev_timem;
00564
00565 mefel_right_hand_side(lcid, fp, fl);
00566 Mp->time = btime;
00567
00568 Gtm->update_elem(Mp->time);
00569
00570 Gtm->update_nodes();
00571
00572
00573 Mt->comreac();
00574
00575 Mt->elemprescdisp();
00576
00577 Mt->clean_ip_new_elem();
00578 }
00579
00580
00581 approximation_temper ();
00582
00583 approximation_humid ();
00584
00585
00586 mefel_right_hand_side (lcid,f,fl);
00587
00588
00589
00590 Mb->comp_sum (f);
00591 Mb->comp_sum_react();
00592 Mb->comp_sum_pdreact();
00593 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fx=%le, fy=%le, fz=%le",
00594 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00595 fflush(Out);
00596
00597 if ((mnce - mnae) > 0){
00598
00599
00600
00601
00602 Mp->strcomp = 0;
00603 internal_forces(lcid,fi);
00604 Mp->strcomp = 1;
00605
00606
00607 subv(fp, fi, fb, nm);
00608 }
00609 else
00610 nullv(fb, nm);
00611
00612
00613 incr_internal_forces (lcid,fi);
00614
00615
00616 for (j=0;j<nm;j++){
00617 fb[j]=f[j]-fp[j]+fi[j];
00618 }
00619
00620
00621 if ((Mp->nlman->stmat==tangent_stiff) || (ii == 1) || (Smat==NULL))
00622 stiffness_matrix (lcid);
00623
00624
00625 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00626
00627
00628 addv(r, dr, nm);
00629
00630
00631 internal_forces (lcid,fi);
00632
00633
00634 norfa=Psolm->pss (f,f,Out);
00635
00636
00637 subv(fl, fi, fb, nm);
00638
00639
00640 norf=Psolm->pss (fb,fb,Out);
00641 if (norfa > 0.0)
00642 norf /= norfa;
00643
00644 if ((Myrank==0)&&(Mespr==1))
00645 fprintf (stdout,"\n\n Norf before inner iteration loop, norf = %e\n\n\n",norf);
00646 j = 0;
00647 if (Psolm->compare_on_master(norf, err) <= 0)
00648 {
00649
00650 nsts++;
00651
00652 nii=1;
00653 }
00654 else
00655 nii=0;
00656
00657 if (nii==0){
00658
00659
00660 fillm(0.0, lsm_a);
00661 fillv(0.0, lsm_r);
00662 for (j=0;j<ini;j++)
00663 {
00664 if((Cp->fcsolv == fullnewtonc) ||
00665 ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0)))
00666 {
00667
00668 stiffness_matrix (lcid);
00669 }
00670
00671 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00672
00673
00674 for (k=0;k<nm;k++)
00675 r[k]+=dr[k];
00676
00677
00678 internal_forces (lcid,fi);
00679
00680
00681 subv(fl, fi, fb, nm);
00682
00683
00684 norf=Psolm->pss (fb,fb,Out);
00685 if (norfa > 0.0)
00686 norf /= norfa;
00687
00688
00689 if ((Myrank==0)&&(Mespr==1))
00690 fprintf (stdout,"\n Inner loop number j = %ld, norf = %e \n",j,norf);
00691 if (Psolm->compare_on_master(norf, err) <= 0)
00692 {
00693
00694
00695 nsts++;
00696
00697 nii=1;
00698 break;
00699 }
00700 else
00701 nii=0;
00702
00703
00704 if (j > 1)
00705 {
00706 lsm_res = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,1);
00707 if (Psolm->compare_on_master(lsm_res, 0.0) > 0)
00708 {
00709 if ((Myrank==0) && (Mespr==1))
00710 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00711
00712 nii=0;
00713 break;
00714 }
00715 }
00716 else
00717 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
00718 }
00719 }
00720
00721 if(Cp->fcsolv == fullnewtonc){
00722
00723 if (Smat != NULL){
00724 delete Smat;
00725 Smat=NULL;
00726 }
00727 }
00728
00729 if (nii==0)
00730 {
00731
00732 if ((Myrank==0) && (Mespr==1))
00733 fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00734 mefel_convergence = 0;
00735
00736
00737 for (j=0;j<nm;j++)
00738 r[j]=lhsb[j];
00739
00740
00741
00742 dt/=2.0;
00743 Mp->timecon.oldtime ();
00744 ii--;
00745
00746 if ((Myrank==0) && (Mespr==1))
00747 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00748 if (Psolm->compare_on_master(dt, dtmin) < 0){
00749 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00750 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
00751 compute_req_val (lcid);
00752 print_step_forced(lcid, i, Mp->time, f);
00753 print_flush();
00754 break;
00755 }
00756 }
00757 else
00758 {
00759
00760 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00761
00762 mefel_convergence = 1;
00763
00764
00765 for (j=0;j<nm;j++)
00766 fp[j]=f[j];
00767
00768
00769 Mt->save_nodval(lcid);
00770
00771
00772 Gtm->update_auxinf();
00773 prev_timem = Mp->time;
00774
00775 dtdef = Mp->timecon.actualforwtimeincr ();
00776 if (nsts==2){
00777 dt*=2.0;
00778 nsts=0;
00779 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00780 }
00781 if (dt>dtdef)
00782 dt=dtdef;
00783
00784
00785 Mb->comp_sum (fi);
00786 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fix=%le, fiy=%le, fiz=%le",
00787 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00788 fflush(Out);
00789
00790
00791 Mm->updateipval();
00792 compute_req_val (lcid);
00793 print_step(lcid, ii, Mp->time, f);
00794 print_flush();
00795
00796 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.hdbtype))
00797 {
00798 if (Mespr==1)
00799 fprintf (stdout,"\n Creating backup file for MEFEL\n");
00800 solver_save (r,f,ii,Mp->time,dt,&Mp->timecon,Ndofm);
00801 }
00802 }
00803 }
00804
00805 if(mefel_convergence == 0){
00806
00807 for (j=0;j<nt;j++){
00808 lhst[j]=lhstb[j];
00809 tdlhst[j]=tdlhstb[j];
00810 }
00811 }
00812 else{
00813 print_stept(0,i,Tp->time,NULL);
00814 print_flusht();
00815 if ((Tp->timecont.isitimptime()==1) && Tp->hdbcont.save_stat())
00816 {
00817 if (Mesprt==1)
00818 fprintf (stdout,"\n Creating backup file for TRFEL\n");
00819 solvert_save (lhst,tdlhst,rhst,i,Cp->time,dt,Cp->timecon,Ndoft);
00820 }
00821 }
00822 } while(Cp->time<=end_time);
00823
00824 print_close ();
00825 print_closet ();
00826
00827 delete [] fi;
00828 delete [] fb;
00829 delete [] fp;
00830 delete [] dr;
00831 delete [] fl;
00832 delete [] lhsb;
00833 delete [] ifn;
00834
00835 delete [] p;
00836 delete [] d;
00837
00838 delete [] lhstb;
00839 delete [] tdlhstb;
00840 }
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852 void par_newton_raphson_gparcoupl_nonlin (long lcid)
00853 {
00854 int mefel_convergence;
00855 long i,ii,j,k,nm,nt,ini,init,nsts,mncd,tncd,mnce,tnce,nii,mnae,nifn;
00856 long *ifn;
00857 double zero,dt,dtmin,dtdef,end_time,alpha,prev_timem,prev_timet,btime;
00858 double *dr,*fb,*r,*d,*f,*fp,*fi,*fl,*p,*lhst,*tdlhst,*rhst,*lhsb,*lhstb,*tdlhstb;
00859 double *fbt,*fit,*lhst_last;
00860
00861 double norf_last, lsm_res;
00862 double norf, norfa,err,errt;
00863 matrix lsm_a(3,3);
00864 vector lsm_r(3), lsm_l(3);
00865
00866
00867 ini = Mp->niilnr;
00868
00869 err = Mp->errnr;
00870
00871 init = Tp->nii;
00872
00873 errt = Tp->err;
00874
00875 nm=Ndofm;
00876
00877 nt=Ndoft;
00878
00879
00880
00881
00882
00883 lhst=Lsrst->give_lhs (lcid);
00884
00885 tdlhst=Lsrst->give_tdlhs (lcid);
00886
00887 rhst=Lsrst->give_rhs (lcid);
00888
00889
00890 d = new double [nt];
00891 nullv (d,nt);
00892
00893 p = new double [nt];
00894 nullv (p,nt);
00895 fbt = new double [nt];
00896 nullv (fbt,nt);
00897 fit = new double [nt];
00898 nullv (fit,nt);
00899 lhst_last = new double [nt];
00900
00901
00902 lhstb = new double [nt];
00903 nullv (lhstb,nt);
00904
00905 tdlhstb = new double [nt];
00906 nullv (tdlhstb,nt);
00907
00908
00909
00910
00911
00912
00913 r = Lsrs->give_lhs (lcid);
00914
00915 f = Lsrs->give_rhs (lcid);
00916
00917
00918 dr = new double [nm];
00919 nullv (dr,nm);
00920
00921 fp = new double [nm];
00922 nullv (fp,nm);
00923
00924 fl = new double [nm];
00925 nullv (fl,nm);
00926
00927 fi = new double [nm];
00928 nullv (fi,nm);
00929
00930 fb = new double [nm];
00931 nullv (fb,nm);
00932
00933 lhsb = new double [nm];
00934 nullv (lhsb,nm);
00935
00936 ifn = new long[Mt->nn];
00937 memset(ifn, 0, sizeof(*ifn)*Mt->nn);
00938
00939
00940 alpha=Tp->alpha;
00941 zero=Tp->zero;
00942
00943
00944
00945 Cp->time=Cp->timecon.starttime ();
00946
00947 dt=Cp->timecon.initialtimeincr ();
00948
00949 end_time = Cp->timecon.endtime ();
00950
00951
00952 Mp->time = Cp->time;
00953 Mp->timecon.take_values (Cp->timecon);
00954
00955 dtmin=Mp->timecon.dtmin;
00956
00957 prev_timem = Mp->time;
00958
00959 Tp->time = Cp->time;
00960 Tp->timecont.take_values (Cp->timecon);
00961
00962 prev_timet = Tp->time;
00963
00964
00965 approximation_inittemper ();
00966
00967
00968 approximation ();
00969
00970 actual_previous_nodval ();
00971
00972
00973
00974
00975 i=0;
00976 ii=0;
00977
00978 nsts=0;
00979
00980
00981 if (Mp->hdbcont.restore_stat()){
00982 if (Mespr==1)
00983 fprintf (stdout,"\n Reading of backup file for MEFEL\n");
00984 solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm);
00985 Mt->rhs_save (fp);
00986
00987 Mm->initmaterialmodels(lcid);
00988 print_init(-1, "at",Pcp->fnim,Pcp->feim);
00989
00990
00991 }
00992 else
00993 {
00994
00995 Mm->initmaterialmodels(lcid);
00996 print_init(-1, "wt",Pcp->fnim,Pcp->feim);
00997 print_step(lcid, i, Mp->time, f);
00998 print_flush();
00999 }
01000 if (Tp->hdbcont.restore_stat()){
01001 if (Mesprt==1)
01002 fprintf (stdout,"\n Reading of backup file for TRFEL\n");
01003 solvert_restore (lhst,tdlhst,f,i,Cp->time,dt,Cp->timecon,Ndoft);
01004
01005 Tm->initmaterialmodels ();
01006 compute_req_valt (0);
01007 print_initt(-1, "at",Pcp->fnit,Pcp->feit);
01008
01009
01010 }
01011 else{
01012
01013 Tm->initmaterialmodels ();
01014 compute_req_valt (0);
01015 print_initt(-1, "wt",Pcp->fnit,Pcp->feit);
01016 print_stept(0,i,Tp->time,NULL);
01017 print_flusht();
01018 }
01019
01020
01021 mefel_convergence = 1;
01022
01023
01024 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01025
01026 Mt->save_nodval(lcid);
01027
01028
01029
01030
01031
01032 do{
01033 if(mefel_convergence == 0){
01034 Cp->timecon.oldtime ();
01035 Tp->timecont.oldtime ();
01036 if (Mesprc==1) fprintf (stdout,"\n\n --------------------------------------------------------------------------------------------");
01037 if (Mesprc==1) fprintf (stdout,"\n REDUCED TIME STEP, METR TIME STEP = %ld : METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
01038 if (Mesprc==1) fprintf (stdout,"\n --------------------------------------------------------------------------------------------");
01039 i--;
01040 }
01041
01042
01043
01044 Cp->time = Cp->timecon.newtime (dt);
01045 Tp->time = Tp->timecont.newtime (dt);
01046 Tp->time = Cp->time;
01047 Tp->timecont.take_values (Cp->timecon);
01048 i++;
01049
01050 if(mefel_convergence == 1){
01051 if (Mesprc==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------------------");
01052 if (Mesprc==1) fprintf (stdout,"\n NEW METR TIME STEP = %ld: METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
01053 if (Mesprc==1) fprintf (stdout,"\n -------------------------------------------------------------------------------------------");
01054 }
01055 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
01056 if (Mesprt==1) fprintf (stdout,"\n TRFEL TIME STEP = %ld, TRFEL TIME = %e, trfel time increment = %e",i,Tp->time,dt);
01057 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
01058
01059
01060
01061
01062 tnce = Gtt->search_newelem (Tp->time,prev_timet);
01063
01064
01065
01066
01067 if (tnce>0){
01068 Tt->initial_nodval ();
01069 }
01070
01071
01072
01073
01074 tncd = Gtt->search_newdofs (Tp->time,prev_timet);
01075
01076
01077 Gtt->update_elem (Tp->time);
01078
01079
01080 Gtt->update_nodes ();
01081
01082
01083 Gtt->update_dofs (Tp->time);
01084
01085
01086
01087 nt=Ndoft;
01088
01089
01090 if (tncd>0){
01091 if (Kmat != NULL){
01092 delete Kmat;
01093 Kmat=NULL;
01094 }
01095
01096 if (Cmat != NULL){
01097 delete Cmat;
01098 Cmat=NULL;
01099 }
01100 }
01101
01102 if (tnce!=0 || tncd!=0){
01103
01104 Tt->lhs_restore (lhst,Lsrst->lhsi,tdlhst);
01105 }
01106
01107 Tm->updateipval ();
01108
01109
01110 conductivity_matrix (lcid);
01111
01112
01113 capacity_matrix (lcid);
01114
01115
01116 for (j=0;j<nt;j++){
01117 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
01118 }
01119
01120
01121 trfel_right_hand_side (0,rhst,nt);
01122
01123
01124 Kmat->gmxv (p,d);
01125
01126
01127
01128 Kmat->scalgm (dt*alpha);
01129 Kmat->addgm (1.0,*Cmat);
01130
01131 for (j=0;j<nt;j++){
01132
01133 rhst[j] = rhst[j] - d[j];
01134
01135
01136 }
01137
01138
01139 Psolt->par_linear_solver (Gtt,Kmat,tdlhst,rhst,Outt,Mesprt);
01140
01141
01142 for (j=0;j<nt;j++){
01143 lhstb[j]=lhst[j];
01144 tdlhstb[j]=tdlhst[j];
01145 }
01146
01147
01148 for (j=0;j<nt;j++){
01149
01150
01151
01152
01153 lhst[j] = p[j] + alpha*dt*tdlhst[j];
01154 }
01155
01156
01157 solution_correction ();
01158
01159 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01160
01161
01162 approximation ();
01163
01164 norf_last = 1.0e20;
01165
01166 for (j=0;j<init;j++){
01167
01168 solution_correction ();
01169
01170 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01171
01172
01173 approximation ();
01174
01175
01176 if (Tp->trsolv == fullnewtont){
01177
01178
01179 capacity_matrix (0);
01180
01181
01182 conductivity_matrix (0);
01183
01184
01185 Kmat->gmxv (p,d);
01186
01187
01188
01189 Kmat->scalgm (dt*alpha);
01190 Kmat->addgm (1.0,*Cmat);
01191 }
01192
01193 if (Tp->trestype==lrhst){
01194
01195 trfel_right_hand_side (0,rhst,nt);
01196
01197 Kmat->gmxv (tdlhst,fit);
01198
01199 for (k=0;k<nt;k++){
01200 fbt[k] = rhst[k] - d[k] - fit[k];
01201 }
01202 }
01203
01204 if (Tp->trestype==fluxest){
01205
01206 internal_fluxes (fit,nt);
01207
01208 for (k=0;k<nt;k++){
01209 fbt[k]=fit[k];
01210 }
01211 }
01212
01213
01214 norf=Psolt->pss (fbt,fbt, Out);
01215
01216 if ((Myrank==0)&&(Mesprt==1))
01217 fprintf (stdout,"\n TRFEL inner loop number j = %ld, norf = %e",j,norf);
01218
01219 if (Psolt->compare_on_master(norf, errt) <= 0)
01220 break;
01221
01222 Psolt->par_linear_solver (Gtt,Kmat,d,fbt,Outt,Mesprt);
01223
01224 for (k=0;k<nt;k++){
01225 tdlhst[k]+=d[k];
01226 lhst[k]+=alpha*dt*d[k];
01227 }
01228
01229
01230 solution_correction ();
01231
01232 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01233
01234
01235 approximation ();
01236
01237
01238
01239 if (Tp->convergcontrolt==yes){
01240 if (Psolt->compare_on_master(norf, norf_last)>0){
01241 for (k=0;k<nt;k++){
01242 lhst[k] = lhst_last[k];
01243 }
01244
01245
01246 solution_correction ();
01247
01248 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01249
01250
01251 approximation ();
01252
01253 if ((Myrank==0)&&(Mesprt==1)) fprintf (stdout,"\n\nTRFEL convergence control: inner iteration skipped %ld error %e\n\n",j,norf);
01254
01255 break;
01256 }
01257 for (k=0;k<nt;k++){
01258 lhst_last[k] = lhst[k];
01259 }
01260 norf_last = norf;
01261 }
01262 }
01263
01264 if (Cp->fcsolv == fullnewtonc){
01265
01266
01267 if (Kmat != NULL){
01268 delete Kmat;
01269 Kmat=NULL;
01270 }
01271
01272 if (Cmat != NULL){
01273 delete Cmat;
01274 Cmat=NULL;
01275 }
01276 }
01277
01278
01279
01280
01281 mefel_convergence = 1;
01282
01283 Mp->time = Mp->timecon.newtime (dt);
01284
01285 if (Mp->time <= Tp->time){
01286
01287 ii++;
01288
01289
01290 Mp->strainstate=0;
01291
01292
01293 Mp->stressstate=0;
01294
01295
01296 Mp->otherstate=0;
01297
01298
01299 for (j=0;j<nm;j++){
01300 lhsb[j]=r[j];
01301 }
01302
01303 if (Mespr==1) fprintf (stdout,"\n\n -----------------------------------------------------------------------------------");
01304 if (Mespr==1) fprintf (stdout,"\n MEFEL TIME STEP = %ld, MEFEL TIME = %e, mefel time increment = %e",ii,Mp->time,dt);
01305 if (Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------------------------");
01306
01307
01308
01309 mnce = Gtm->search_changed_elem(Mp->time, mnae);
01310
01311
01312 mncd = Gtm->search_changed_dofs(Mp->time,prev_timem);
01313
01314 nifn = 0;
01315 if (mnce!=0){
01316
01317
01318 nifn = Gtm->search_iface_nodes(ifn);
01319 }
01320
01321 if ((mnce!=0) || (mncd!=0))
01322 {
01323
01324
01325 Mt->save_node_inidispl(lcid, Mp->time, prev_timem);
01326
01327
01328
01329
01330
01331 if (nifn)
01332 {
01333
01334 Gtm->switch_new_elem();
01335
01336
01337 Gtm->update_nodes();
01338
01339
01340
01341 Gtm->update_active_dofs(Mp->time);
01342
01343
01344 Ndofm = Gtm->codenum_generation(Out);
01345
01346
01347
01348 Mt->restore_nodval(r, Ndofm);
01349
01350
01351 stress_initdispl(lcid);
01352
01353
01354
01355 Gtm->update_active_dofs(Mp->time);
01356
01357 Gtm->clear_intf_dofs(ifn);
01358
01359
01360 Ndofm = Gtm->codenum_generation(Out);
01361
01362
01363 Mp->strcomp = 0;
01364 internal_forces(lcid,fb);
01365 cmulv(-1.0, fb, Ndofm);
01366 Mp->strcomp = 1;
01367
01368
01369
01370 delete Smat;
01371 Smat=NULL;
01372 stiffness_matrix (lcid);
01373 nullv(r, Ndofm);
01374
01375
01376 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
01377
01378
01379 addv(r, dr, nm);
01380
01381
01382
01383 Mt->save_elem_inidispl(lcid, ifn);
01384
01385
01386
01387 Gtm->remove_nodes(ifn);
01388
01389
01390
01391 Mt->save_nodval(lcid);
01392 }
01393
01394
01395 Gtm->update_elem(Mp->time);
01396
01397
01398 Gtm->update_nodes();
01399
01400
01401 Gtm->update_dofs (Mp->time);
01402
01403
01404 nm = Ndofm = Gtm->codenum_generation(Out);
01405
01406
01407 delete Smat;
01408 Smat=NULL;
01409
01410
01411 Mt->restore_nodval(r,nm);
01412
01413
01414
01415
01416
01417
01418 Gtm->update_elem(prev_timem);
01419
01420 Gtm->update_nodes();
01421
01422 btime = Mp->time;
01423
01424 Mp->time = prev_timem;
01425
01426 mefel_right_hand_side(lcid, fp, fl);
01427 Mp->time = btime;
01428
01429 Gtm->update_elem(Mp->time);
01430
01431 Gtm->update_nodes();
01432
01433
01434 Mt->comreac();
01435
01436 Mt->elemprescdisp();
01437
01438 Mt->clean_ip_new_elem();
01439 }
01440
01441
01442 approximation_temper ();
01443
01444 approximation_humid ();
01445
01446
01447 mefel_right_hand_side (lcid,f,fl);
01448
01449
01450
01451 Mb->comp_sum (f);
01452 Mb->comp_sum_react();
01453 Mb->comp_sum_pdreact();
01454 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fx=%le, fy=%le, fz=%le",
01455 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
01456 fflush(Out);
01457
01458 if ((mnce - mnae) > 0)
01459 {
01460
01461
01462
01463
01464 Mp->strcomp = 0;
01465 internal_forces(lcid,fi);
01466 Mp->strcomp = 1;
01467
01468
01469 subv(fp, fi, fb, nm);
01470 }
01471 else
01472 nullv(fb, nm);
01473
01474
01475
01476 incr_internal_forces (lcid,fi);
01477
01478
01479 for (j=0;j<nm;j++){
01480 fb[j]=f[j]-fp[j]+fi[j];
01481 }
01482
01483
01484 if ((Mp->nlman->stmat==tangent_stiff) || (ii == 1) || (Smat==NULL))
01485 stiffness_matrix (lcid);
01486
01487
01488 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
01489
01490
01491 addv(r, dr, nm);
01492
01493
01494 internal_forces (lcid,fi);
01495
01496
01497 norfa=Psolm->pss (f,f,Out);
01498
01499
01500 subv(fl, fi, fb, nm);
01501
01502
01503 norf=Psolm->pss (fb,fb,Out);
01504 if (norfa > 0.0)
01505 norf /= norfa;
01506
01507 if ((Myrank==0)&&(Mespr==1)) fprintf (stdout,"\n\n Norf before inner iteration loop norf = %e\n\n\n",norf);
01508 j = 0;
01509 if (Psolm->compare_on_master(norf, err) <= 0)
01510 {
01511
01512 nsts++;
01513
01514 nii=1;
01515 }
01516 else
01517 nii=0;
01518
01519 if (nii==0){
01520
01521
01522 fillm(0.0, lsm_a);
01523 fillv(0.0, lsm_r);
01524 for (j=0;j<ini;j++)
01525 {
01526 if((Cp->fcsolv == fullnewtonc) ||
01527 ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0)))
01528 {
01529
01530 stiffness_matrix (lcid);
01531 }
01532
01533 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
01534
01535
01536 for (k=0;k<nm;k++)
01537 r[k]+=dr[k];
01538
01539
01540 internal_forces (lcid,fi);
01541
01542
01543 subv(fl, fi, fb, nm);
01544
01545
01546 norf=Psolm->pss (fb,fb,Out);
01547 if (norfa > 0.0)
01548 norf /= norfa;
01549
01550
01551 if ((Myrank==0)&&(Mespr==1))
01552 fprintf (stdout,"\n Inner loop j = %ld norf = %e\n",j,norf);
01553
01554 if (Psolm->compare_on_master(norf, err) <= 0)
01555 {
01556
01557
01558 nsts++;
01559
01560 nii=1;
01561 break;
01562 }
01563 else
01564 nii=0;
01565
01566
01567 if (j > 1)
01568 {
01569 lsm_res = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,1);
01570 if (Psolm->compare_on_master(lsm_res, 0.0) > 0)
01571 {
01572 if ((Myrank==0)&&(Mespr==1))
01573 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
01574
01575 nii=0;
01576 break;
01577 }
01578 }
01579 else
01580 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
01581 }
01582 }
01583
01584 if(Cp->fcsolv == fullnewtonc){
01585
01586 if (Smat != NULL){
01587 delete Smat;
01588 Smat=NULL;
01589 }
01590 }
01591
01592 if (nii==0)
01593 {
01594
01595 if ((Myrank==0) && (Mespr==1))
01596 fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
01597
01598 mefel_convergence = 0;
01599
01600
01601 for (j=0;j<nm;j++)
01602 r[j]=lhsb[j];
01603
01604
01605
01606 dt/=2.0;
01607 Mp->timecon.oldtime ();
01608 ii--;
01609
01610 if ((Myrank==0) && (Mespr==1))
01611 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
01612 if (Psolm->compare_on_master(dt, dtmin) < 0){
01613 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n\n time increment is less than minimum time increment");
01614 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
01615 compute_req_val (lcid);
01616 print_step_forced(lcid, i, Mp->time, f);
01617 print_flush();
01618 break;
01619 }
01620 }
01621 else
01622 {
01623
01624 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
01625 mefel_convergence = 1;
01626
01627
01628 for (j=0;j<nm;j++)
01629 fp[j]=f[j];
01630
01631
01632 Mt->save_nodval(lcid);
01633 prev_timem = Mp->time;
01634
01635 dtdef = Mp->timecon.actualforwtimeincr ();
01636 if (nsts==2){
01637 dt*=2.0;
01638 nsts=0;
01639 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
01640 }
01641 if (dt>dtdef)
01642 dt=dtdef;
01643
01644
01645 Mb->comp_sum (fi);
01646 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fix=%le, fiy=%le, fiz=%le",
01647 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
01648 fflush(Out);
01649
01650
01651 Mm->updateipval();
01652 compute_req_val (lcid);
01653 print_step(lcid, ii, Mp->time, f);
01654 print_flush();
01655
01656 if ((Mp->timecon.isitimptime ()==1) && (Mp->hdbcont.hdbtype))
01657 {
01658 if (Mespr==1)
01659 fprintf (stdout,"\n Creating backup file for MEFEL\n");
01660 solver_save (r,f,ii,Mp->time,dt,&Mp->timecon,Ndofm);
01661 }
01662 }
01663 }
01664 if(mefel_convergence == 0){
01665
01666 for (j=0;j<nt;j++){
01667 lhst[j]=lhstb[j];
01668 tdlhst[j]=tdlhstb[j];
01669 }
01670 }
01671 else{
01672 print_stept(0,i,Tp->time,NULL);
01673 print_flusht();
01674 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat()){
01675 if (Mesprt==1)
01676 fprintf (stdout,"\n Creating backup file for TRFEL\n");
01677 solvert_save (lhst,tdlhst,rhst,i,Cp->time,dt,Cp->timecon,Ndoft);
01678 }
01679 }
01680 }while(Cp->time<=end_time);
01681
01682 print_close ();
01683 print_closet ();
01684
01685 delete [] fi;
01686 delete [] fb;
01687 delete [] fp;
01688 delete [] dr;
01689 delete [] fl;
01690 delete [] lhsb;
01691 delete [] ifn;
01692
01693 delete [] p;
01694 delete [] d;
01695
01696 delete [] fit;
01697 delete [] fbt;
01698 delete [] lhst_last;
01699
01700 delete [] lhstb;
01701 delete [] tdlhstb;
01702 }