00001 #include <math.h>
00002 #include <string.h>
00003 #include "ppcsolver.h"
00004 #include "pglobal.h"
00005 #include "pglobalt.h"
00006 #include "pglobalc.h"
00007 #include "globalc.h"
00008 #include "global.h"
00009 #include "globalt.h"
00010 #include "globmat.h"
00011 #include "globmatt.h"
00012 #include "globmatc.h"
00013 #include "elemswitcht.h"
00014 #include "gmatrix.h"
00015 #include "gtopology.h"
00016 #include "mechprint.h"
00017 #include "transprint.h"
00018 #include "backupsol.h"
00019 #include "backupsolt.h"
00020 #include "seqfilesm.h"
00021 #include "seqfilest.h"
00022 #include "genfile.h"
00023 #include "mpi.h"
00024
00025
00026
00027
00028
00029
00030
00031
00032 void par_solve_pcouplprob ()
00033 {
00034 long lcid;
00035
00036
00037
00038 lcid=0;
00039
00040 switch (Cp->tnlinsol){
00041 case newtonc:{
00042
00043 par_newton_raphson_parcoupl_comp (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_parcoupl_comp (long lcid)
00063 {
00064 int mefel_convergence,trfel_convergence;
00065 double dt,end_time;
00066 long i,k,nt,ti,tli,tnsts,tret;
00067 double tnewtime,tdt, tdtmin,tdtmax, tdtdef, tend_time;
00068 double *lhsb,*tdlhsb;
00069 np_glob_vec np_gv;
00070 long tlcid = 0;
00071 long mi, mli, mnsts, mret;
00072 double mnewtime,mdt, mdtmin,mdtmax, mdtdef, mend_time;
00073 mt_glob_vec mt_gv;
00074 long mlcid = 0;
00075
00076
00077
00078
00079
00080 Cp->time = Cp->timecon.starttime ();
00081
00082 dt = Cp->timecon.initialtimeincr ();
00083
00084 end_time = Cp->timecon.endtime ();
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 init_trfel_mefel ();
00098
00099
00100
00101
00102
00103 par_nonstat_trfel_init(tlcid, np_gv);
00104
00105 tdt = Tp->timecont.initialtimeincr ();
00106
00107 tdtmin=Tp->timecont.dtmin;
00108
00109 tdtmax=Tp->timecont.dtmax;
00110
00111 tend_time = Tp->timecont.endtime ();
00112
00113 tli = ti = np_gv.istep;
00114
00115 tret = 0;
00116
00117
00118
00119
00120
00121 par_visco_mefel_init(mlcid, mt_gv);
00122
00123 mdt = Mp->timecon.initialtimeincr ();
00124
00125 mdtmin=Mp->timecon.dtmin;
00126
00127 mdtmax=Mp->timecon.dtmax;
00128
00129 mend_time = Mp->timecon.endtime ();
00130
00131
00132 mli = mi = mt_gv.istep;
00133
00134 mnsts=0;
00135
00136 mret = 0;
00137
00138
00139 if(Tp->timecont.tct > 0){
00140 nt=Ndoft;
00141 lhsb = new double [nt];
00142 tdlhsb = new double [nt];
00143 nullv (lhsb,nt);
00144 nullv (tdlhsb,nt);
00145 }
00146
00147
00148
00149
00150 i=0;
00151 mefel_convergence = 1;
00152 trfel_convergence = 1;
00153
00154 do{
00155 i++;
00156 if(mefel_convergence == 0 || trfel_convergence == 0){
00157 Tp->timecont.oldtime ();
00158
00159 if (Myrank==0 && Mesprc==1)
00160 fprintf (stdout,"\n\n --------------------------------------------------------------------------------------------");
00161 if (Myrank==0 && Mesprc==1)
00162 fprintf (stdout,"\n REDUCED TIME STEP, PMETR TIME STEP = %ld : PMETR TIME = %e, pmetr time increment = %e",i,Cp->time,dt);
00163 if (Myrank==0 && Mesprc==1)
00164 fprintf (stdout,"\n --------------------------------------------------------------------------------------------");
00165 i--;
00166 }
00167 else{
00168
00169
00170
00171 tnewtime = Tp->time = Tp->timecont.newtime (tdt);
00172 Cp->time = Tp->time;
00173
00174 if (Myrank==0 && Mesprc==1)
00175 fprintf (stdout,"\n\n -------------------------------------------------------------------------------------");
00176 if (Myrank==0 && Mesprc==1)
00177 fprintf (stdout,"\n NEW PMETR TIME STEP = %ld: PMETR TIME = %e, pmetr time increment = %e",i,Cp->time,dt);
00178 if (Myrank==0 && Mesprc==1)
00179 fprintf (stdout,"\n -------------------------------------------------------------------------------------");
00180 }
00181
00182
00183 copy_data();
00184
00185
00186 trfel_convergence = 1;
00187 ti = i;
00188
00189
00190 if((Tp->timecont.timefun.tfunc != stat) && Tp->timecont.tct == 0)
00191 tdt = Tp->timecont.actualbacktimeincr ();
00192
00193 if(Tp->tprob == nonstationary_problem)
00194
00195
00196
00197 tret = par_one_step_trfel_linear(tlcid, tnewtime, tdt, ti, tli, np_gv);
00198
00199 if(Tp->tprob == nonlinear_nonstationary_problem)
00200
00201
00202
00203 tret = par_one_step_trfel_nonlinear(tlcid, tnewtime, tdt, ti, tli, np_gv);
00204
00205
00206 if(Tp->tprob == nonlinear_nonstationary_problem && Tp->timecont.tct > 0){
00207 if (tret >= 0)
00208 {
00209 trfel_convergence = 1;
00210
00211 for (k=0;k<nt;k++){
00212 lhsb[k] = np_gv.lhs[k];
00213 tdlhsb[k] = np_gv.tdlhs[k];
00214 }
00215
00216 if (tret == 0)
00217 tnsts++;
00218 tdtdef = Tp->timecont.actualforwtimeincr();
00219 if (tnsts==2)
00220 {
00221 tdt*=2.0;
00222 tnsts=0;
00223
00224 if (Mesprt==1)
00225 fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00226 }
00227 if (tdt<=tdtdef)
00228 tdt = tdtdef;
00229
00230 if (tdt>tdtmax)
00231 tdt = tdtmax;
00232 }
00233 else
00234 {
00235
00236
00237 trfel_convergence = 0;
00238
00239 tdt/=2.0;
00240
00241 if (Mesprt==1)
00242 fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00243
00244 if (tdt<tdtmin)
00245 {
00246 if (Mesprt==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00247 if (Mesprt==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00248 if (Mesprt==1) fprintf (stderr,"\n FORCED output of results from this step is performed\n");
00249 compute_req_valt (lcid);
00250 print_stept(lcid, ti, Tp->time, np_gv.rhs);
00251 print_flusht();
00252 break;
00253 }
00254 }
00255 }
00256
00257
00258 if (Cp->cleanmatrix == clean_yes){
00259
00260
00261 if (Kmat != NULL){
00262 fprintf (stderr,"\n\n Cleaning of conductivity matrix \n\n ");
00263 delete Kmat;
00264 Kmat=NULL;
00265 }
00266
00267 if (Cmat != NULL){
00268 fprintf (stderr,"\n\n Cleaning of capacity matrix \n\n ");
00269 delete Cmat;
00270 Cmat=NULL;
00271 }
00272 }
00273
00274
00275
00276 mefel_convergence = 1;
00277
00278 if((Mp->timecon.timefun.tfunc != stat) && Mp->timecon.tct == 0)
00279 mdt = Mp->timecon.actualbacktimeincr ();
00280
00281
00282
00283 if (Mp->time < Tp->time)
00284 mnewtime = Mp->time = Mp->timecon.newtime (mdt);
00285
00286 if (Mp->time <= Tp->time){
00287
00288 mi++;
00289
00290
00291 copy_data();
00292
00293
00294
00295 mret = par_one_step_mefel(mlcid, mnewtime, mdt, mi, mli, mt_gv);
00296
00297
00298 if(Mp->timecon.tct > 0){
00299 if (mret >= 0)
00300 {
00301 mefel_convergence = 1;
00302
00303 if (mret == 0)
00304 mnsts++;
00305 mdtdef = Mp->timecon.actualforwtimeincr();
00306 if (mnsts==2)
00307 {
00308 mdt*=2.0;
00309 mnsts=0;
00310 if (Myrank==0 && Mespr==1)
00311 fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00312 }
00313 if (mdt>mdtdef)
00314 mdt=mdtdef;
00315
00316 if (mdt>mdtmax)
00317 mdt = mdtmax;
00318 }
00319 else
00320 {
00321
00322
00323 mefel_convergence = 0;
00324
00325 mdt/=2.0;
00326 Mp->timecon.oldtime ();
00327 mi--;
00328
00329 if (Myrank==0 && Mespr==1)
00330 fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00331 if (mdt<mdtmin)
00332 {
00333 if (Myrank==0 && Mespr==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00334 if (Myrank==0 && Mespr==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00335 if (Myrank==0 && Mespr==1) fprintf (stderr,"\n FORCED output of results from this step is performed\n");
00336 compute_req_val (mlcid);
00337 print_step_forced(mlcid, mi, Mp->time, mt_gv.fb);
00338
00339 print_flush();
00340 break;
00341 }
00342 }
00343 }
00344
00345
00346 if(Cp->cleanmatrix == clean_yes){
00347
00348 if (Smat != NULL){
00349 fprintf (stderr,"\n\n Cleaning of stiffness matrix \n\n ");
00350 delete Smat;
00351 Smat=NULL;
00352 }
00353 }
00354 }
00355
00356 }while(Cp->time<=end_time);
00357 print_close ();
00358 print_closet();
00359
00360 if(Tp->timecont.tct > 0){
00361 delete [] lhsb;
00362 delete [] tdlhsb;
00363 }
00364 }
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374 void par_nonstat_trfel_init (long lcid, np_glob_vec &np_gv)
00375 {
00376 long n;
00377 double dt;
00378
00379
00380 n=Ndoft;
00381
00382 np_gv.alloc(n);
00383
00384 np_gv.lhs = Lsrst->give_lhs (lcid);
00385
00386 np_gv.tdlhs = Lsrst->give_tdlhs (lcid);
00387
00388 np_gv.rhs = Lsrst->give_rhs (lcid);
00389
00390 if(Tp->tprob == nonlinear_nonstationary_problem)
00391 np_gv.alloc_aux(n);
00392
00393
00394 Tp->time=Tp->timecont.starttime ();
00395
00396 dt = Tp->timecont.initialtimeincr ();
00397
00398
00399 approximation();
00400
00401
00402 if (Tp->nvs==1 && Tp->pnvs==1)
00403 actual_previous_nodval ();
00404
00405
00406 if (Tp->hdbcont.restore_stat()){
00407 solvert_restore (np_gv.lhs, np_gv.tdlhs, np_gv.f, np_gv.istep, Tp->time, dt, Tp->timecont, n);
00408
00409 Tm->initmaterialmodels();
00410 compute_req_valt (lcid);
00411 print_initt(-1, "at",Pcp->fnit,Pcp->feit);
00412 }
00413 else{
00414
00415 Tm->initmaterialmodels();
00416 compute_req_valt (lcid);
00417 print_initt(-1, "wt",Pcp->fnit,Pcp->feit);
00418 print_stept(lcid, np_gv.istep, Tp->time, np_gv.rhs);
00419 }
00420 print_flusht();
00421
00422 }
00423
00424
00425
00426
00427
00428
00429
00430 void par_visco_mefel_init(long lcid, mt_glob_vec &mt_gv)
00431 {
00432 long n;
00433 double dt;
00434
00435 n=Ndofm;
00436
00437 mt_gv.alloc(n);
00438
00439 mt_gv.r = Lsrs->give_lhs (lcid);
00440
00441 mt_gv.f = Lsrs->give_rhs (lcid);
00442
00443
00444 Mp->time=Mp->timecon.starttime ();
00445
00446
00447 if (Mp->hdbcont.restore_stat())
00448 {
00449 if (Mespr==1)
00450 fprintf (stdout,"\n Reading of backup file\n");
00451 solver_restore (mt_gv.r, mt_gv.fp, mt_gv.istep, Mp->time, dt, &Mp->timecon, n);
00452
00453 Mm->initmaterialmodels(lcid);
00454 print_init(-1, "at",Pcp->fnim,Pcp->feim);
00455
00456 print_flush();
00457 }
00458 else
00459 {
00460 mt_gv.istep=0;
00461
00462 Mm->initmaterialmodels(lcid);
00463
00464 if (Mp->eigstrains > 0)
00465 internal_forces(lcid, mt_gv.fp);
00466
00467 print_init(-1, "wt",Pcp->fnim,Pcp->feim);
00468 print_step(lcid, mt_gv.istep, Mp->time, mt_gv.f);
00469
00470 print_flush();
00471 }
00472 }
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483 long par_one_step_trfel_linear (long lcid,double time, double dt, long istep, long li, np_glob_vec &np_gv)
00484 {
00485 long j,n;
00486 double alpha,*f,*d,*p,*lhs,*tdlhs,*rhs;
00487
00488
00489 Tp->time = time;
00490
00491 Tp->istep = istep;
00492
00493
00494 lhs = np_gv.lhs;
00495
00496 tdlhs = np_gv.tdlhs;
00497
00498 rhs = np_gv.rhs;
00499
00500 f = np_gv.f;
00501
00502 d = np_gv.d;
00503
00504 p = np_gv.p;
00505
00506
00507 alpha=Tp->alpha;
00508
00509
00510 n=Ndoft;
00511
00512 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n\n ------------------------------------------------------------------------");
00513 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n PTRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
00514 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------\n");
00515
00516
00517 Tm->updateipval ();
00518
00519
00520 capacity_matrix (lcid);
00521
00522
00523 conductivity_matrix (lcid);
00524
00525 if (istep == 0){
00526 Psolt->computation_stat (Gtt,Kmat);
00527 }
00528
00529
00530
00531 for (j=0;j<n;j++){
00532 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00533 }
00534
00535
00536
00537 Kmat->gmxv (d,p);
00538
00539
00540
00541 Kmat->scalgm (dt*alpha);
00542 Kmat->addgm (1.0,*Cmat);
00543
00544
00545 trfel_right_hand_side (0,rhs,n);
00546
00547
00548
00549 for (j=0;j<n;j++){
00550 f[j] = rhs[j] - p[j];
00551 }
00552
00553
00554
00555
00556 Psolt->par_linear_solver (Gtt,Kmat,tdlhs,f,Outt,Mesprt);
00557
00558
00559
00560 for (j=0;j<n;j++){
00561 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00562 }
00563
00564 solution_correction ();
00565
00566 approximation ();
00567
00568
00569 if (Tp->nvs==1 && Tp->pnvs==1)
00570 actual_previous_nodval ();
00571
00572 compute_req_valt (0);
00573 print_stept(0,istep,Tp->time,rhs);
00574 print_flusht();
00575
00576 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00577 solvert_save (lhs,tdlhs,f,istep,Tp->time,dt,Tp->timecont,n);
00578
00579 if (j == 0)
00580 return 0;
00581
00582 return 1;
00583 }
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594 long par_one_step_trfel_nonlinear (long lcid,double time, double dt, long istep, long li, np_glob_vec &np_gv)
00595 {
00596 long j,k,n,ini;
00597 double alpha,*f,*d,*p,*lhs,*tdlhs,*rhs;
00598 double *fb,*fi;
00599 double norf_last;
00600 double zero,norf,*err,*thresh;
00601
00602
00603 Tp->time = time;
00604
00605 Tp->istep = istep;
00606
00607
00608 n=Ndoft;
00609
00610
00611 lhs = np_gv.lhs;
00612
00613 tdlhs = np_gv.tdlhs;
00614
00615 rhs = np_gv.rhs;
00616
00617 f = np_gv.f;
00618
00619 d = np_gv.d;
00620
00621 p = np_gv.p;
00622
00623
00624 fb = np_gv.fb;
00625 fi = np_gv.fi;
00626
00627
00628 nullv (fb,n);
00629 nullv (fi,n);
00630
00631
00632
00633 ini = Tp->nii;
00634
00635 err = Tp->errarr;
00636
00637 thresh=Tp->threshrhs;
00638
00639
00640 alpha=Tp->alpha;
00641
00642 zero=Tp->zero;
00643
00644 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n\n -------------------------------------------------------------------");
00645 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n PTRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
00646 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n -------------------------------------------------------------------\n");
00647
00648
00649 n=Ndoft;
00650
00651
00652 Tm->updateipval ();
00653
00654
00655 capacity_matrix (lcid);
00656
00657
00658 conductivity_matrix (lcid);
00659
00660 if (istep == 0){
00661 Psolt->computation_stat (Gtt,Kmat);
00662 }
00663
00664
00665
00666 for (j=0;j<n;j++){
00667 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00668 }
00669
00670
00671
00672 Kmat->gmxv (d,p);
00673
00674
00675
00676 Kmat->scalgm (dt*alpha);
00677 Kmat->addgm (1.0,*Cmat);
00678 Kmat->copygm (*Cmat);
00679
00680
00681 trfel_right_hand_side (0,rhs,n);
00682
00683
00684
00685 for (j=0;j<n;j++){
00686 f[j] = rhs[j] - p[j];
00687 }
00688
00689
00690
00691
00692 Psolt->par_linear_solver (Gtt,Kmat,tdlhs,f,Outt,Mesprt);
00693
00694
00695
00696 for (j=0;j<n;j++){
00697 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00698 }
00699
00700
00701 solution_correction ();
00702
00703 approximation ();
00704
00705
00706 if (Tp->nvs==1 && Tp->pnvs==1)
00707 actual_previous_nodval ();
00708
00709
00710 norf_last = 1.0e20;
00711
00712 for (j=0;j<ini;j++){
00713
00714
00715 if (Tp->trsolv == fullnewtont){
00716
00717
00718 capacity_matrix (0);
00719
00720
00721 conductivity_matrix (0);
00722
00723
00724 Kmat->gmxv (d,p);
00725
00726
00727
00728 Kmat->scalgm (dt*alpha);
00729 Kmat->addgm (1.0,*Cmat);
00730 Kmat->copygm (*Cmat);
00731 }
00732
00733 if (Tp->trestype==lrhst){
00734
00735 trfel_right_hand_side (lcid,rhs,n);
00736
00737 Kmat->gmxv (tdlhs,fi);
00738
00739 for (k=0;k<n;k++){
00740 fb[k] = rhs[k] - p[k] - fi[k];
00741 }
00742 }
00743
00744 if (Tp->trestype==fluxest){
00745
00746 internal_fluxes (fi,n);
00747
00748 for (k=0;k<n;k++){
00749 fb[k]=fi[k];
00750 }
00751 }
00752
00753 norf = Psolt->unbalanced_values (fb,Outt);
00754
00755 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n inner iteration %ld error %e",j,norf);
00756 if (norf<err[0]) break;
00757
00758 Psolt->par_linear_solver (Gtt,Cmat,fi,fb,Outt,Mesprt);
00759
00760 for (k=0;k<n;k++){
00761 tdlhs[k]+=fi[k];
00762 lhs[k]+=alpha*dt*fi[k];
00763 }
00764
00765
00766 solution_correction ();
00767
00768 approximation ();
00769
00770
00771
00772 if (Tp->nvs==1 && Tp->pnvs==1)
00773 actual_previous_nodval ();
00774
00775
00776
00777 if (Tp->convergcontrolt==yes){
00778 if (norf > norf_last){
00779
00780 if (Myrank==0 && Mesprt==1)
00781 fprintf (stdout,"\n\n convergence control: inner iteration skiped %ld error %e\n\n",j,norf);
00782
00783 return -2;
00784 }
00785
00786 norf_last = norf;
00787 }
00788 }
00789
00790 if(norf >= err[0])
00791 return -1;
00792
00793
00794
00795 compute_req_valt (lcid);
00796 print_stept(lcid,istep,Tp->time,rhs);
00797 print_flusht();
00798
00799 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00800 solvert_save (lhs,tdlhs,f,istep,Tp->time,dt,Tp->timecont,n);
00801
00802 return 0;
00803 }
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828 long par_one_step_mefel (long lcid,double time, double dt, long istep, long li, mt_glob_vec &mt_gv)
00829 {
00830 long j,n,ini;
00831 double *f,*fl,*fi,*fb,*fp,*r,*dr,*lhsb;
00832 double norf, err, zero;
00833 matrix lsm_a(3,3);
00834 vector lsm_r(3), lsm_l(3);
00835 Mp->nlman->divc_step = new double[10];
00836 Mp->nlman->divc_err = new double[10];
00837
00838
00839 Mp->time = time;
00840
00841 Mp->istep = istep;
00842
00843
00844 Mp->strainstate=0;
00845
00846
00847 Mp->stressstate=0;
00848
00849
00850 Mp->otherstate=0;
00851
00852
00853 r = mt_gv.r;
00854
00855 f = mt_gv.f;
00856
00857 dr = mt_gv.dr;
00858
00859 fp = mt_gv.fp;
00860
00861 fl = mt_gv.fl;
00862
00863 fi = mt_gv.fi;
00864
00865 fb = mt_gv.fb;
00866
00867 lhsb = mt_gv.lhsb;
00868
00869
00870 n=Ndofm;
00871
00872
00873 ini = Mp->niilnr;
00874
00875
00876 err = Mp->errnr;
00877
00878
00879 zero = Mp->zero;
00880
00881
00882 for (j=0;j<n;j++)
00883 lhsb[j]=r[j];
00884
00885
00886 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n\n ------------------------------------------------------------------------");
00887 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n PMEFEL Time step = %ld, Time %e, Time increment = %e",istep,Mp->time,dt);
00888 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n ------------------------------------------------------------------------\n");
00889
00890
00891 mefel_right_hand_side (lcid,f,fl);
00892
00893
00894
00895 Mb->comp_sum (f);
00896
00897
00898 fflush(Out);
00899
00900
00901 incr_internal_forces (lcid,fi);
00902
00903
00904
00905 subv(f, fp, fb, n);
00906 addv(fb, fi, fb, n);
00907
00908 norf = par_gnewton_raphson_one_step_mefel(lcid, Mp->nlman, fl, r, fb, dr, fi, istep-1, j, li, ini, err);
00909
00910 if (norf>err)
00911 {
00912
00913 if (Myrank==0 && Mespr==1)
00914 fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00915
00916
00917
00918 copyv(lhsb, r, n);
00919
00920 return -1;
00921 }
00922
00923
00924 if ((Myrank==0 && Mespr==1) && (j > 0))
00925 fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00926
00927
00928
00929 copyv(f, fp, n);
00930
00931
00932 Mb->comp_sum (fi);
00933
00934
00935 fflush(Out);
00936
00937
00938 Mm->updateipval();
00939 compute_req_val (lcid);
00940 print_step(lcid, istep, Mp->time, fl);
00941
00942 print_flush();
00943
00944 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.save_stat()))
00945 {
00946 if (Mespr==1)
00947 fprintf (stdout,"\n Creating backup file\n");
00948 solver_save (r,fp,istep,Mp->time,dt,&Mp->timecon,n);
00949 }
00950 if (j == 0)
00951 return 0;
00952
00953 return 1;
00954 }
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977 double par_gnewton_raphson_one_step_mefel(long lcid, nonlinman *nlman, double *fa, double *ra, double *fb, double *dr, double *fi,
00978 long istep, long &j, long li, long ini, double ierr)
00979 {
00980 long n = Ndofm;
00981 double norf, norfa;
00982 matrix lsm_a(3,3);
00983 vector lsm_r(3), lsm_l(3);
00984
00985
00986 assemble_stiffness_matrix(lcid,istep,-1,li);
00987
00988
00989 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00990
00991
00992
00993 addv(ra, dr, n);
00994
00995
00996 internal_forces (lcid,fi);
00997
00998
00999
01000 subv(fa, fi, fb, n);
01001
01002
01003 norfa=Psolm->pss (fa,fa,Out);
01004
01005 norf=Psolm->pss(fb,fb,Out);
01006 if (norfa != 0.0)
01007 norf /= norfa;
01008
01009 if ((Myrank==0)&&(Mespr==1))
01010 fprintf (stdout,", norm %e",norf);
01011
01012 if (norf<ierr)
01013 {
01014 j = -1;
01015 return norf;
01016 }
01017
01018
01019 fillm(0.0, lsm_a);
01020 fillv(0.0, lsm_r);
01021 for (j=0; j<ini; j++)
01022 {
01023 Mp->jstep = j;
01024
01025
01026 assemble_stiffness_matrix(lcid,istep,j,li);
01027
01028
01029 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
01030
01031
01032 addv(ra, dr, n);
01033
01034
01035 internal_forces(lcid,fi);
01036
01037
01038
01039 subv(fa, fi, fb, n);
01040
01041
01042 norf=Psolm->pss(fb,fb,Out);
01043 if (norfa != 0.0)
01044 norf /= norfa;
01045
01046 if ((Myrank==0) && (Mespr==1))
01047 fprintf (stdout,"\n increment %ld inner loop j %ld norf %e",istep,j,norf);
01048
01049 if (norf < ierr)
01050
01051 return norf;
01052
01053
01054 if (check_divergency(nlman, lsm_a, lsm_r, lsm_l, j, norf))
01055 return norf;
01056 }
01057
01058
01059 return norf;
01060 }
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075 void par_newton_raphson_parcoupl (long lcid)
01076 {
01077 switch (Tp->tprob){
01078 case nonstationary_problem:{
01079
01080
01081 par_newton_raphson_parcoupl_lin (lcid);
01082 break;
01083 }
01084 case nonlinear_nonstationary_problem:{
01085
01086
01087 par_newton_raphson_parcoupl_nonlin (lcid);
01088 break;
01089 }
01090 default:{
01091 print_err("unknown METR problem type is required", __FILE__,__LINE__,__func__);
01092 }
01093 }
01094 }
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106 void par_newton_raphson_parcoupl_lin (long lcid)
01107 {
01108 int mefel_convergence;
01109 long i,ii,j,k,nm,nt,ini,nsts,nii;
01110 double zero,dt,dtmin,dtdef,end_time,alpha,s;
01111 double *dr,*fb,*r,*d,*f,*fp,*fi,*fl,*p,*lhst,*tdlhst,*rhst,*lhsb,*lhstb,*tdlhstb;
01112 double norf, norfa, err, lsm_res;
01113 matrix lsm_a(3,3);
01114 vector lsm_r(3), lsm_l(3);
01115 double alphadtn;
01116 double alphadto = 0.0;
01117 double zerodt = 1.0e-12;
01118
01119
01120 ini = Mp->niilnr;
01121
01122 err = Mp->errnr;
01123
01124 nm=Ndofm;
01125
01126 nt=Ndoft;
01127
01128
01129
01130
01131
01132
01133 lhst=Lsrst->give_lhs (lcid);
01134
01135 tdlhst=Lsrst->give_tdlhs (lcid);
01136
01137 rhst=Lsrst->give_rhs (lcid);
01138
01139
01140 d = new double [nt];
01141 nullv (d,nt);
01142
01143 p = new double [nt];
01144 nullv (p,nt);
01145
01146 lhstb = new double [nt];
01147 nullv (lhstb,nt);
01148
01149 tdlhstb = new double [nt];
01150 nullv (tdlhstb,nt);
01151
01152
01153
01154
01155
01156
01157 r = Lsrs->give_lhs (lcid);
01158
01159 f = Lsrs->give_rhs (lcid);
01160
01161
01162 dr = new double [nm];
01163 nullv (dr,nm);
01164
01165 fp = new double [nm];
01166 nullv (fp,nm);
01167
01168 fl = new double [nm];
01169 nullv (fl,nm);
01170
01171 fi = new double [nm];
01172 nullv (fi,nm);
01173
01174 fb = new double [nm];
01175 nullv (fb,nm);
01176
01177 lhsb = new double [nm];
01178 nullv (lhsb,nm);
01179
01180
01181 alpha=Tp->alpha;
01182 zero=Tp->zero;
01183
01184
01185
01186 Cp->time=Cp->timecon.starttime ();
01187
01188 dt=Cp->timecon.initialtimeincr ();
01189
01190 end_time = Cp->timecon.endtime ();
01191
01192
01193 Mp->time = Cp->time;
01194 Mp->timecon.take_values (Cp->timecon);
01195
01196 dtmin=Mp->timecon.dtmin;
01197
01198 Tp->time = Cp->time;
01199 Tp->timecont.take_values (Cp->timecon);
01200
01201
01202 approximation_inittemper();
01203
01204 approximation();
01205
01206 actual_previous_nodval();
01207
01208
01209
01210
01211 i=0;
01212 ii=0;
01213
01214 nsts=0;
01215
01216
01217 if (Mp->hdbcont.restore_stat()){
01218 if (Mespr==1)
01219 fprintf (stdout,"\n Reading of backup file for MEFEL\n");
01220 solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm);
01221
01222 Mm->initmaterialmodels(lcid);
01223 print_init(-1, "at",Pcp->fnim,Pcp->feim);
01224
01225
01226 }
01227 else
01228 {
01229
01230 Mm->initmaterialmodels(lcid);
01231 print_init(-1, "wt",Pcp->fnim,Pcp->feim);
01232 print_step(lcid, i, Mp->time, f);
01233 print_flush();
01234 }
01235 if (Tp->hdbcont.restore_stat()){
01236 if (Mesprt==1)
01237 fprintf (stdout,"\n Reading of backup file for TRFEL\n");
01238 solvert_restore (lhst,tdlhst,f,i,Cp->time,dt,Cp->timecon,Ndoft);
01239
01240 Tm->initmaterialmodels();
01241 compute_req_valt (0);
01242 print_initt(-1, "at",Pcp->fnit,Pcp->feit);
01243
01244
01245 }
01246 else{
01247
01248 Tm->initmaterialmodels();
01249 compute_req_valt (0);
01250 print_initt(-1, "wt",Pcp->fnit,Pcp->feit);
01251 print_stept(0,i,Tp->time,NULL);
01252 print_flusht();
01253 }
01254
01255 mefel_convergence = 1;
01256
01257
01258
01259
01260 do{
01261 if(mefel_convergence == 0){
01262 Cp->timecon.oldtime ();
01263 Tp->timecont.oldtime ();
01264 if ((Myrank==0) && (Mesprc==1)) fprintf (stdout,"\n\n --------------------------------------------------------------------------------------------");
01265 if ((Myrank==0) && (Mesprc==1)) fprintf (stdout,"\n REDUCED TIME STEP, METR TIME STEP = %ld : METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
01266 if ((Myrank==0) && (Mesprc==1)) fprintf (stdout,"\n --------------------------------------------------------------------------------------------");
01267 i--;
01268 }
01269
01270
01271
01272 Cp->time = Cp->timecon.newtime (dt);
01273 Tp->time = Tp->timecont.newtime (dt);
01274 Tp->time = Cp->time;
01275 Tp->timecont.take_values (Cp->timecon);
01276 i++;
01277
01278 if(mefel_convergence == 1){
01279 if ((Myrank==0) && (Mesprc==1)) fprintf (stdout,"\n\n -------------------------------------------------------------------------------------------");
01280 if ((Myrank==0) && (Mesprc==1)) fprintf (stdout,"\n NEW METR TIME STEP = %ld: METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
01281 if ((Myrank==0) && (Mesprc==1)) fprintf (stdout,"\n -------------------------------------------------------------------------------------------");
01282 }
01283 if ((Myrank==0) && (Mesprt==1)) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
01284 if ((Myrank==0) && (Mesprt==1)) fprintf (stdout,"\n TRFEL TIME STEP = %ld, TRFEL TIME = %e, trfel time increment = %e",i,Tp->time,dt);
01285 if ((Myrank==0) && (Mesprt==1)) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
01286
01287 Tm->updateipval ();
01288
01289 alphadtn = alpha*dt;
01290 if (Myrank==0)
01291 {
01292 fprintf(stdout, "\n adto=%le, adtn=%le, adto/adtn-1=%le, Cmat=%p, Kmat=%p, fcsolv=%d\n",
01293 alphadto, alphadtn, fabs(alphadto/alphadtn - 1.0), Cmat, Kmat, Cp->fcsolv);
01294 fflush(stdout);
01295 }
01296 if ((fabs(alphadto/alphadtn - 1.0) > zerodt) || (fabs(alphadtn) < zero) ||
01297 (Cmat==NULL) || (Kmat == NULL) || (Cp->fcsolv != linearc))
01298 {
01299
01300 conductivity_matrix (lcid);
01301
01302
01303 capacity_matrix (lcid);
01304 }
01305
01306 for (j=0;j<nt;j++){
01307 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
01308 }
01309
01310
01311 trfel_right_hand_side (0,rhst,nt);
01312
01313
01314 Kmat->gmxv (p,d);
01315
01316
01317 if ((fabs(alphadto/alphadtn - 1.0) > zerodt) || (fabs(alphadtn) < zero) ||
01318 (Cmat==NULL) || (Kmat == NULL) || (Cp->fcsolv != linearc))
01319 {
01320
01321
01322
01323
01324
01325 Cmat->addgm(dt*alpha, *Kmat);
01326 }
01327
01328 for (j=0;j<nt;j++){
01329 rhst[j] = rhst[j] - d[j];
01330 d[j]=tdlhst[j];
01331 }
01332
01333
01334
01335 Psolt->par_linear_solver (Gtt,Cmat,tdlhst,rhst,Outt,Mesprt);
01336
01337
01338 for (j=0;j<nt;j++){
01339 lhstb[j]=lhst[j];
01340 tdlhstb[j]=tdlhst[j];
01341 }
01342
01343
01344 for (j=0;j<nt;j++){
01345 s=(1.0-alpha)*d[j]+alpha*tdlhst[j];
01346 lhst[j]+=dt*s;
01347 }
01348
01349
01350 solution_correction ();
01351
01352 approximation ();
01353
01354 actual_previous_nodval ();
01355
01356 if (Cp->fcsolv == fullnewtonc){
01357
01358
01359 if (Kmat != NULL){
01360 delete Kmat;
01361 Kmat=NULL;
01362 }
01363
01364 if (Cmat != NULL){
01365 delete Cmat;
01366 Cmat=NULL;
01367 }
01368 }
01369
01370
01371
01372 mefel_convergence = 1;
01373
01374 Mp->time = Mp->timecon.newtime (dt);
01375
01376 if (Mp->time <= Tp->time){
01377
01378 ii++;
01379
01380
01381 Mp->strainstate=0;
01382
01383
01384 Mp->stressstate=0;
01385
01386
01387 Mp->otherstate=0;
01388
01389
01390 approximation_temper ();
01391
01392 approximation_humid ();
01393
01394
01395 for (j=0;j<nm;j++){
01396 lhsb[j]=r[j];
01397 }
01398
01399 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n -------------------------------------------------------------------------------");
01400 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n MEFEL TIME STEP = %ld, MEFEL TIME = %e, mefel time increment = %e",ii,Mp->time,dt);
01401 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n -------------------------------------------------------------------------------");
01402
01403
01404 mefel_right_hand_side (lcid,f,fl);
01405
01406
01407
01408 Mb->comp_sum (f);
01409
01410
01411 fflush(Out);
01412
01413
01414 incr_internal_forces (lcid,fi);
01415
01416
01417 for (j=0;j<nm;j++){
01418 fb[j]=f[j]-fp[j]+fi[j];
01419 }
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431 if ((Mp->nlman->stmat==tangent_stiff) || (ii == 1) || (Smat==NULL))
01432 stiffness_matrix (lcid);
01433
01434
01435 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
01436
01437
01438 for (j=0;j<nm;j++){
01439 r[j]+=dr[j];
01440 }
01441
01442
01443 internal_forces (lcid,fi);
01444
01445
01446 norfa=Psolm->pss (f,f,Out);
01447
01448 for (j=0;j<nm;j++){
01449 fb[j] = fl[j] - fi[j];
01450 }
01451
01452 norf=Psolm->pss (fb,fb,Out);
01453 if (norfa > 0.0)
01454 norf /= norfa;
01455
01456 if ((Myrank==0)&&(Mespr==1))
01457 fprintf (stdout,"\n\n Norf before inner iteration loop norf = %e\n\n\n",norf);
01458 j = 0;
01459 if (Psolm->compare_on_master(norf, err) <= 0)
01460 {
01461
01462 nsts++;
01463
01464 nii=1;
01465 }
01466 else
01467 nii=0;
01468
01469 if (nii==0){
01470
01471
01472 fillm(0.0, lsm_a);
01473 fillv(0.0, lsm_r);
01474 for (j=0;j<ini;j++)
01475 {
01476 if((Cp->fcsolv == fullnewtonc) ||
01477 ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0)))
01478 {
01479
01480 stiffness_matrix (lcid);
01481 }
01482
01483 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
01484
01485
01486 for (k=0;k<nm;k++)
01487 r[k]+=dr[k];
01488
01489
01490 internal_forces (lcid,fi);
01491
01492
01493 for (k=0;k<nm;k++)
01494 fb[k]= fl[k] - fi[k];
01495
01496
01497 norf=Psolm->pss (fb,fb,Out);
01498 if (norfa > 0.0)
01499 norf /= norfa;
01500
01501
01502 if ((Myrank==0)&&(Mespr==1))
01503 fprintf (stdout,"\n Inner loop number j = %ld, norf = %e \n",j,norf);
01504 if (Psolm->compare_on_master(norf, err) <= 0)
01505 {
01506
01507
01508 nsts++;
01509
01510 nii=1;
01511 break;
01512 }
01513 else
01514 nii=0;
01515
01516
01517 if (j > 1)
01518 {
01519 lsm_res = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,1);
01520 if (Psolm->compare_on_master(lsm_res, 0.0) > 0)
01521 {
01522 if ((Myrank==0) && (Mespr==1))
01523 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
01524
01525 nii=0;
01526 break;
01527 }
01528 }
01529 else
01530 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
01531 }
01532 }
01533
01534 if(Cp->fcsolv == fullnewtonc){
01535
01536 if (Smat != NULL){
01537 delete Smat;
01538 Smat=NULL;
01539 }
01540 }
01541
01542 if (nii==0)
01543 {
01544
01545 if ((Myrank==0) && (Mespr==1))
01546 fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
01547 mefel_convergence = 0;
01548
01549
01550 for (j=0;j<nm;j++)
01551 r[j]=lhsb[j];
01552
01553
01554
01555 dt/=2.0;
01556 Mp->timecon.oldtime ();
01557 ii--;
01558
01559 if ((Myrank==0) && (Mespr==1))
01560 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
01561 if (Psolm->compare_on_master(dt, dtmin) < 0){
01562 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n\n time increment is less than minimum time increment");
01563 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
01564 compute_req_val (lcid);
01565 print_step_forced(lcid, i, Mp->time, f);
01566 print_flush();
01567 break;
01568 }
01569 }
01570 else
01571 {
01572
01573 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
01574 mefel_convergence = 1;
01575
01576
01577 for (j=0;j<nm;j++)
01578 fp[j]=f[j];
01579
01580 dtdef = Mp->timecon.actualforwtimeincr ();
01581 if (nsts==2){
01582 dt*=2.0;
01583 nsts=0;
01584 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
01585 }
01586 if (dt>dtdef)
01587 dt=dtdef;
01588
01589
01590 Mb->comp_sum (fi);
01591
01592
01593 fflush(Out);
01594
01595
01596 Mm->updateipval();
01597 compute_req_val (lcid);
01598 print_step(lcid, ii, Mp->time, f);
01599 print_flush();
01600
01601 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.hdbtype))
01602 {
01603 if (Mespr==1)
01604 fprintf (stdout,"\n Creating backup file for MEFEL\n");
01605 solver_save (r,f,ii,Mp->time,dt,&Mp->timecon,Ndofm);
01606 }
01607 }
01608 }
01609
01610 if(mefel_convergence == 0){
01611
01612 for (j=0;j<nt;j++){
01613 lhst[j]=lhstb[j];
01614 tdlhst[j]=tdlhstb[j];
01615 }
01616 }
01617 else{
01618 alphadto = alphadtn;
01619 print_stept(0,i,Tp->time,NULL);
01620 print_flusht();
01621 if ((Tp->timecont.isitimptime()==1) && Tp->hdbcont.save_stat())
01622 {
01623 if (Mesprt==1)
01624 fprintf (stdout,"\n Creating backup file for TRFEL\n");
01625 solvert_save (lhst,tdlhst,rhst,i,Cp->time,dt,Cp->timecon,Ndoft);
01626 }
01627 }
01628 } while(Cp->time<=end_time);
01629
01630 print_close ();
01631 print_closet ();
01632
01633 delete [] fi;
01634 delete [] fb;
01635 delete [] fp;
01636 delete [] dr;
01637 delete [] fl;
01638 delete [] lhsb;
01639
01640 delete [] p;
01641 delete [] d;
01642
01643 delete [] lhstb;
01644 delete [] tdlhstb;
01645 }
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657 void par_newton_raphson_parcoupl_nonlin (long lcid)
01658 {
01659 int mefel_convergence;
01660 long i,ii,j,k,nm,nt,ini,init,nsts,nii;
01661 double zero,dt,dtmin,dtdef,end_time,alpha;
01662 double *dr,*fb,*r,*d,*f,*fp,*fi,*fl,*p,*lhst,*tdlhst,*rhst,*lhsb,*lhstb,*tdlhstb;
01663 double *fbt,*fit,*lhst_last;
01664
01665 double norf_last, lsm_res;
01666 double norf, norfa,err,errt;
01667 matrix lsm_a(3,3);
01668 vector lsm_r(3), lsm_l(3);
01669
01670
01671 ini = Mp->niilnr;
01672
01673 err = Mp->errnr;
01674
01675 init = Tp->nii;
01676
01677 errt = Tp->err;
01678
01679 nm=Ndofm;
01680
01681 nt=Ndoft;
01682
01683
01684
01685
01686
01687 lhst=Lsrst->give_lhs (lcid);
01688
01689 tdlhst=Lsrst->give_tdlhs (lcid);
01690
01691 rhst=Lsrst->give_rhs (lcid);
01692
01693
01694 d = new double [nt];
01695 nullv (d,nt);
01696
01697 p = new double [nt];
01698 nullv (p,nt);
01699 fbt = new double [nt];
01700 nullv (fbt,nt);
01701 fit = new double [nt];
01702 nullv (fit,nt);
01703 lhst_last = new double [nt];
01704
01705
01706 lhstb = new double [nt];
01707 nullv (lhstb,nt);
01708
01709 tdlhstb = new double [nt];
01710 nullv (tdlhstb,nt);
01711
01712
01713
01714
01715
01716
01717 r = Lsrs->give_lhs (lcid);
01718
01719 f = Lsrs->give_rhs (lcid);
01720
01721
01722 dr = new double [nm];
01723 nullv (dr,nm);
01724
01725 fp = new double [nm];
01726 nullv (fp,nm);
01727
01728 fl = new double [nm];
01729 nullv (fl,nm);
01730
01731 fi = new double [nm];
01732 nullv (fi,nm);
01733
01734 fb = new double [nm];
01735 nullv (fb,nm);
01736
01737 lhsb = new double [nm];
01738 nullv (lhsb,nm);
01739
01740
01741 alpha=Tp->alpha;
01742 zero=Tp->zero;
01743
01744
01745
01746 Cp->time=Cp->timecon.starttime ();
01747
01748 dt=Cp->timecon.initialtimeincr ();
01749
01750 end_time = Cp->timecon.endtime ();
01751
01752
01753 Mp->time = Cp->time;
01754 Mp->timecon.take_values (Cp->timecon);
01755
01756 dtmin=Mp->timecon.dtmin;
01757
01758 Tp->time = Cp->time;
01759 Tp->timecont.take_values (Cp->timecon);
01760
01761
01762 approximation_inittemper ();
01763
01764
01765 approximation ();
01766
01767 actual_previous_nodval ();
01768
01769
01770
01771
01772 i=0;
01773 ii=0;
01774
01775 nsts=0;
01776
01777
01778 if (Mp->hdbcont.restore_stat()){
01779 if (Mespr==1)
01780 fprintf (stdout,"\n Reading of backup file for MEFEL\n");
01781 solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm);
01782
01783 Mm->initmaterialmodels(lcid);
01784 print_init(-1, "at",Pcp->fnim,Pcp->feim);
01785
01786
01787 }
01788 else
01789 {
01790
01791 Mm->initmaterialmodels(lcid);
01792 print_init(-1, "wt",Pcp->fnim,Pcp->feim);
01793 print_step(lcid, i, Mp->time, f);
01794 print_flush();
01795 }
01796 if (Tp->hdbcont.restore_stat()){
01797 if (Mesprt==1)
01798 fprintf (stdout,"\n Reading of backup file for TRFEL\n");
01799 solvert_restore (lhst,tdlhst,f,i,Cp->time,dt,Cp->timecon,Ndoft);
01800
01801 Tm->initmaterialmodels();
01802 compute_req_valt (0);
01803 print_initt(-1, "at",Pcp->fnit,Pcp->feit);
01804
01805
01806 }
01807 else{
01808
01809 Tm->initmaterialmodels();
01810 compute_req_valt (0);
01811 print_initt(-1, "wt",Pcp->fnit,Pcp->feit);
01812 print_stept(0,i,Tp->time,NULL);
01813 print_flusht();
01814 }
01815
01816 mefel_convergence = 1;
01817
01818
01819
01820
01821 do{
01822 if(mefel_convergence == 0){
01823 Cp->timecon.oldtime ();
01824 Tp->timecont.oldtime ();
01825 if (Mesprc==1) fprintf (stdout,"\n\n --------------------------------------------------------------------------------------------");
01826 if (Mesprc==1) fprintf (stdout,"\n REDUCED TIME STEP, METR TIME STEP = %ld : METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
01827 if (Mesprc==1) fprintf (stdout,"\n --------------------------------------------------------------------------------------------");
01828 i--;
01829 }
01830
01831
01832
01833 Cp->time = Cp->timecon.newtime (dt);
01834 Tp->time = Tp->timecont.newtime (dt);
01835 Tp->time = Cp->time;
01836 Tp->timecont.take_values (Cp->timecon);
01837 i++;
01838
01839 if(mefel_convergence == 1){
01840 if (Mesprc==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------------------");
01841 if (Mesprc==1) fprintf (stdout,"\n NEW METR TIME STEP = %ld: METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
01842 if (Mesprc==1) fprintf (stdout,"\n -------------------------------------------------------------------------------------------");
01843 }
01844 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
01845 if (Mesprt==1) fprintf (stdout,"\n TRFEL TIME STEP = %ld, TRFEL TIME = %e, trfel time increment = %e",i,Tp->time,dt);
01846 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
01847
01848 Tm->updateipval ();
01849
01850
01851 conductivity_matrix (lcid);
01852
01853
01854 capacity_matrix (lcid);
01855
01856
01857 for (j=0;j<nt;j++){
01858 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
01859 }
01860
01861
01862 trfel_right_hand_side (0,rhst,nt);
01863
01864
01865 Kmat->gmxv (p,d);
01866
01867
01868
01869 Kmat->scalgm (dt*alpha);
01870 Kmat->addgm (1.0,*Cmat);
01871
01872 for (j=0;j<nt;j++){
01873
01874 rhst[j] = rhst[j] - d[j];
01875
01876
01877 }
01878
01879
01880 Psolt->par_linear_solver (Gtt,Kmat,tdlhst,rhst,Outt,Mesprt);
01881
01882
01883 for (j=0;j<nt;j++){
01884 lhstb[j]=lhst[j];
01885 tdlhstb[j]=tdlhst[j];
01886 }
01887
01888
01889 for (j=0;j<nt;j++){
01890
01891
01892
01893
01894 lhst[j] = p[j] + alpha*dt*tdlhst[j];
01895 }
01896
01897
01898 solution_correction ();
01899
01900 approximation ();
01901
01902 norf_last = 1.0e20;
01903
01904 for (j=0;j<init;j++){
01905
01906 solution_correction ();
01907
01908 approximation ();
01909
01910
01911 if (Tp->trsolv == fullnewtont){
01912
01913
01914 capacity_matrix (0);
01915
01916
01917 conductivity_matrix (0);
01918
01919
01920 Kmat->gmxv (p,d);
01921
01922
01923
01924 Kmat->scalgm (dt*alpha);
01925 Kmat->addgm (1.0,*Cmat);
01926 }
01927
01928 if (Tp->trestype==lrhst){
01929
01930 trfel_right_hand_side (0,rhst,nt);
01931
01932 Kmat->gmxv (tdlhst,fit);
01933
01934 for (k=0;k<nt;k++){
01935 fbt[k] = rhst[k] - d[k] - fit[k];
01936 }
01937 }
01938
01939 if (Tp->trestype==fluxest){
01940
01941 internal_fluxes (fit,nt);
01942
01943 for (k=0;k<nt;k++){
01944 fbt[k]=fit[k];
01945 }
01946 }
01947
01948
01949 norf=Psolt->pss (fbt,fbt, Out);
01950
01951 if ((Myrank==0)&&(Mesprt==1))
01952 fprintf (stdout,"\n TRFEL inner loop number j = %ld, norf = %e",j,norf);
01953
01954 if (Psolt->compare_on_master(norf, errt) <= 0)
01955 break;
01956
01957 Psolt->par_linear_solver (Gtt,Kmat,d,fbt,Outt,Mesprt);
01958
01959 for (k=0;k<nt;k++){
01960 tdlhst[k]+=d[k];
01961 lhst[k]+=alpha*dt*d[k];
01962 }
01963
01964
01965 solution_correction ();
01966
01967 approximation ();
01968
01969
01970
01971 if (Tp->convergcontrolt==yes){
01972 if (Psolt->compare_on_master(norf, norf_last)>0){
01973 for (k=0;k<nt;k++){
01974 lhst[k] = lhst_last[k];
01975 }
01976
01977
01978 solution_correction ();
01979
01980 approximation ();
01981
01982 if ((Myrank==0)&&(Mesprt==1)) fprintf (stdout,"\n\nTRFEL convergence control: inner iteration skipped %ld error %e\n\n",j,norf);
01983
01984 break;
01985 }
01986 for (k=0;k<nt;k++){
01987 lhst_last[k] = lhst[k];
01988 }
01989 norf_last = norf;
01990 }
01991 }
01992
01993 if (Cp->fcsolv == fullnewtonc){
01994
01995
01996 if (Kmat != NULL){
01997 delete Kmat;
01998 Kmat=NULL;
01999 }
02000
02001 if (Cmat != NULL){
02002 delete Cmat;
02003 Cmat=NULL;
02004 }
02005 }
02006
02007
02008
02009
02010 mefel_convergence = 1;
02011
02012 Mp->time = Mp->timecon.newtime (dt);
02013
02014 if (Mp->time <= Tp->time){
02015
02016 ii++;
02017
02018
02019 Mp->strainstate=0;
02020
02021
02022 Mp->stressstate=0;
02023
02024
02025 Mp->otherstate=0;
02026
02027
02028 approximation_temper ();
02029
02030 approximation_humid ();
02031
02032
02033 for (j=0;j<nm;j++){
02034 lhsb[j]=r[j];
02035 }
02036
02037 if (Mespr==1) fprintf (stdout,"\n\n -----------------------------------------------------------------------------------");
02038 if (Mespr==1) fprintf (stdout,"\n MEFEL TIME STEP = %ld, MEFEL TIME = %e, mefel time increment = %e",ii,Mp->time,dt);
02039 if (Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------------------------");
02040
02041
02042 mefel_right_hand_side (lcid,f,fl);
02043
02044
02045
02046 Mb->comp_sum (f);
02047 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fx=%le, fy=%le, fz=%le",
02048 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
02049 fflush(Out);
02050
02051
02052 incr_internal_forces (lcid,fi);
02053
02054
02055 for (j=0;j<nm;j++){
02056 fb[j]=f[j]-fp[j]+fi[j];
02057 }
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069 if ((Mp->nlman->stmat==tangent_stiff) || (ii == 1) || (Smat==NULL))
02070 stiffness_matrix (lcid);
02071
02072
02073 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
02074
02075
02076 for (j=0;j<nm;j++){
02077 r[j]+=dr[j];
02078 }
02079
02080
02081 internal_forces (lcid,fi);
02082
02083
02084 norfa=Psolm->pss (f,f,Out);
02085
02086 for (j=0;j<nm;j++){
02087 fb[j] = fl[j] - fi[j];
02088 }
02089
02090 norf=Psolm->pss (fb,fb,Out);
02091 if (norfa > 0.0)
02092 norf /= norfa;
02093
02094 if ((Myrank==0)&&(Mespr==1)) fprintf (stdout,"\n\n Norf before inner iteration loop norf = %e\n\n\n",norf);
02095 j = 0;
02096 if (Psolm->compare_on_master(norf, err) <= 0)
02097 {
02098
02099 nsts++;
02100
02101 nii=1;
02102 }
02103 else
02104 nii=0;
02105
02106 if (nii==0){
02107
02108
02109 fillm(0.0, lsm_a);
02110 fillv(0.0, lsm_r);
02111 for (j=0;j<ini;j++)
02112 {
02113 if((Cp->fcsolv == fullnewtonc) ||
02114 ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0)))
02115 {
02116
02117 stiffness_matrix (lcid);
02118 }
02119
02120 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
02121
02122
02123 for (k=0;k<nm;k++)
02124 r[k]+=dr[k];
02125
02126
02127 internal_forces (lcid,fi);
02128
02129
02130 for (k=0;k<nm;k++)
02131 fb[k]= fl[k] - fi[k];
02132
02133
02134 norf=Psolm->pss (fb,fb,Out);
02135 if (norfa > 0.0)
02136 norf /= norfa;
02137
02138
02139 if ((Myrank==0)&&(Mespr==1))
02140 fprintf (stdout,"\n Inner loop j = %ld norf = %e\n",j,norf);
02141
02142 if (Psolm->compare_on_master(norf, err) <= 0)
02143 {
02144
02145
02146 nsts++;
02147
02148 nii=1;
02149 break;
02150 }
02151 else
02152 nii=0;
02153
02154
02155 if (j > 1)
02156 {
02157 lsm_res = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,1);
02158 if (Psolm->compare_on_master(lsm_res, 0.0) > 0)
02159 {
02160 if (Myrank==0)
02161 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
02162
02163 nii=0;
02164 break;
02165 }
02166 }
02167 else
02168 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
02169 }
02170 }
02171
02172 if(Cp->fcsolv == fullnewtonc){
02173
02174 if (Smat != NULL){
02175 delete Smat;
02176 Smat=NULL;
02177 }
02178 }
02179
02180 if (nii==0)
02181 {
02182
02183 if ((Myrank==0) && (Mespr==1))
02184 fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
02185
02186 mefel_convergence = 0;
02187
02188
02189 for (j=0;j<nm;j++)
02190 r[j]=lhsb[j];
02191
02192
02193
02194 dt/=2.0;
02195 Mp->timecon.oldtime ();
02196 ii--;
02197
02198 if ((Myrank==0) && (Mespr==1))
02199 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
02200 if (Psolm->compare_on_master(dt, dtmin) < 0){
02201 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n\n time increment is less than minimum time increment");
02202 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
02203 compute_req_val (lcid);
02204 print_step_forced(lcid, i, Mp->time, f);
02205 print_flush();
02206 break;
02207 }
02208 }
02209 else
02210 {
02211
02212 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
02213 mefel_convergence = 1;
02214
02215
02216 for (j=0;j<nm;j++)
02217 fp[j]=f[j];
02218
02219 dtdef = Mp->timecon.actualforwtimeincr ();
02220 if (nsts==2){
02221 dt*=2.0;
02222 nsts=0;
02223 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
02224 }
02225 if (dt>dtdef)
02226 dt=dtdef;
02227
02228
02229 Mb->comp_sum (fi);
02230 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fix=%le, fiy=%le, fiz=%le",
02231 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
02232 fflush(Out);
02233
02234
02235 Mm->updateipval();
02236 compute_req_val (lcid);
02237 print_step(lcid, ii, Mp->time, f);
02238 print_flush();
02239
02240 if ((Mp->timecon.isitimptime ()==1) && (Mp->hdbcont.hdbtype))
02241 {
02242 if (Mespr==1)
02243 fprintf (stdout,"\n Creating backup file for MEFEL\n");
02244 solver_save (r,f,ii,Mp->time,dt,&Mp->timecon,Ndofm);
02245 }
02246 }
02247 }
02248 if(mefel_convergence == 0){
02249
02250 for (j=0;j<nt;j++){
02251 lhst[j]=lhstb[j];
02252 tdlhst[j]=tdlhstb[j];
02253 }
02254 }
02255 else{
02256 print_stept(0,i,Tp->time,NULL);
02257 print_flusht();
02258 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat()){
02259 if (Mesprt==1)
02260 fprintf (stdout,"\n Creating backup file for TRFEL\n");
02261 solvert_save (lhst,tdlhst,rhst,i,Cp->time,dt,Cp->timecon,Ndoft);
02262 }
02263 }
02264 }while(Cp->time<=end_time);
02265
02266 print_close ();
02267 print_closet ();
02268
02269 delete [] fi;
02270 delete [] fb;
02271 delete [] fp;
02272 delete [] dr;
02273 delete [] fl;
02274 delete [] lhsb;
02275
02276 delete [] p;
02277 delete [] d;
02278
02279 delete [] fit;
02280 delete [] fbt;
02281 delete [] lhst_last;
02282
02283 delete [] lhstb;
02284 delete [] tdlhstb;
02285 }