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