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