00001 #include "npglvec.h"
00002 #include "backupsolt.h"
00003 #include "transprint.h"
00004 #include "pnpsolvert.h"
00005 #include "pglobalt.h"
00006 #include "seqfilest.h"
00007 #include <string.h>
00008 #include <math.h>
00009 #include "mpi.h"
00010
00011
00012
00013
00014
00015
00016
00017
00018 void par_solve_nonstationary_problem ()
00019 {
00020 if (Tp->tprob == discont_nonstat_problem){
00021
00022
00023 }
00024 else{
00025
00026
00027 par_nonstat_solv_vform_comp ();
00028
00029
00030
00031 }
00032 }
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 void par_nonstat_solv_vform_comp ()
00045 {
00046 long i,k,li,nsts,ret,n;
00047 double newtime,dt,dtmin,dtmax,dtdef,end_time;
00048 double *lhsb,*tdlhsb;
00049 np_glob_vec np_gv;
00050 long lcid = 0;
00051
00052
00053
00054
00055 par_nonstat_solver_init(lcid, np_gv);
00056
00057
00058 dt = Tp->timecont.initialtimeincr ();
00059
00060 dtmin=Tp->timecont.dtmin;
00061
00062 dtmax=Tp->timecont.dtmax;
00063
00064 end_time = Tp->timecont.endtime ();
00065
00066
00067 li = i = np_gv.istep;
00068
00069 nsts=0;
00070
00071 ret = 0;
00072
00073
00074 if(Tp->timecont.tct > 0){
00075 n=Ndoft;
00076 lhsb = new double [n];
00077 tdlhsb = new double [n];
00078 nullv (lhsb,n);
00079 nullv (tdlhsb,n);
00080 }
00081
00082
00083
00084
00085 do{
00086 i++;
00087
00088 newtime = Tp->time = Tp->timecont.newtime();
00089
00090 if((Tp->timecont.timefun.tfunc != stat) && Tp->timecont.tct == 0)
00091 dt = Tp->timecont.actualbacktimeincr ();
00092
00093 if(Tp->tprob == nonstationary_problem)
00094
00095
00096 ret = par_one_step_linear(lcid, newtime, dt, i, li, np_gv);
00097
00098
00099 if(Tp->tprob == nonlinear_nonstationary_problem){
00100
00101
00102 ret = par_one_step_nonlinear(lcid, newtime, dt, i, li, np_gv);
00103 }
00104
00105
00106 if(Tp->tprob == nonlinear_nonstationary_problem && Tp->timecont.tct > 0){
00107 if (ret >= 0){
00108
00109 for (k=0;k<n;k++){
00110 lhsb[k] = np_gv.lhs[k];
00111 tdlhsb[k] = np_gv.tdlhs[k];
00112 }
00113
00114 if (ret == 0)
00115 nsts++;
00116 dtdef = Tp->timecont.actualforwtimeincr();
00117 if (nsts==2)
00118 {
00119 dt*=2.0;
00120 nsts=0;
00121 if (Myrank==0 && Mesprt==1)
00122 fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00123 }
00124 if (dt<=dtdef)
00125 dt = dtdef;
00126
00127 if (dt>dtmax)
00128 dt = dtmax;
00129 }
00130 else{
00131
00132
00133 dt/=2.0;
00134 Tp->timecont.oldtime ();
00135 i--;
00136
00137 for (k=0;k<n;k++){
00138 np_gv.lhs[k] = lhsb[k];
00139 np_gv.tdlhs[k] = tdlhsb[k];
00140 }
00141
00142 if (Myrank==0 && Mesprt==1)
00143 fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00144
00145 if (dt<dtmin)
00146 {
00147 if (Myrank==0 && Mesprt==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00148 if (Myrank==0 && Mesprt==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00149 if (Myrank==0 && Mesprt==1) fprintf (stderr,"\n FORCED output of results from this step is performed\n");
00150 compute_req_valt (lcid);
00151 print_stept(lcid, i, Tp->time, np_gv.rhs);
00152 print_flusht();
00153 break;
00154 }
00155 }
00156 }
00157
00158 }while(Tp->time < end_time);
00159
00160 print_closet();
00161
00162
00163
00164 if(Tp->timecont.tct > 0){
00165 delete [] lhsb;
00166 delete [] tdlhsb;
00167 }
00168 }
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 void par_nonstat_solver_init (long lcid, np_glob_vec &np_gv)
00180 {
00181 long n;
00182 double dt;
00183
00184
00185 n=Ndoft;
00186
00187 np_gv.alloc(n);
00188
00189 np_gv.lhs = Lsrst->give_lhs (lcid);
00190
00191 np_gv.tdlhs = Lsrst->give_tdlhs (lcid);
00192
00193 np_gv.rhs = Lsrst->give_rhs (lcid);
00194
00195 if(Tp->tprob == nonlinear_nonstationary_problem)
00196 np_gv.alloc_aux(n);
00197
00198
00199 Tp->time=Tp->timecont.starttime ();
00200
00201 dt = Tp->timecont.initialtimeincr ();
00202
00203
00204 approximation();
00205
00206
00207 if (Tp->nvs==1 && Tp->pnvs==1)
00208 actual_previous_nodval ();
00209
00210
00211 if (Tp->hdbcont.restore_stat()){
00212 solvert_restore (np_gv.lhs, np_gv.tdlhs, np_gv.f, np_gv.istep, Tp->time, dt, Tp->timecont, n);
00213
00214 Tm->initmaterialmodels();
00215 compute_req_valt (lcid);
00216 print_initt(-1, "at",Ptp->fni,Ptp->fei);
00217 }
00218 else{
00219
00220 Tm->initmaterialmodels();
00221 compute_req_valt (lcid);
00222 print_initt(-1, "wt",Ptp->fni,Ptp->fei);
00223 print_stept(lcid, np_gv.istep, Tp->time, np_gv.rhs);
00224 }
00225 print_flusht();
00226 }
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237 long par_one_step_linear (long lcid,double time, double dt, long istep, long li, np_glob_vec &np_gv)
00238 {
00239 long j,n;
00240 double alpha,*f,*d,*p,*lhs,*tdlhs,*rhs;
00241
00242
00243 Tp->time = time;
00244
00245 Tp->istep = istep;
00246
00247
00248 lhs = np_gv.lhs;
00249
00250 tdlhs = np_gv.tdlhs;
00251
00252 rhs = np_gv.rhs;
00253
00254 f = np_gv.f;
00255
00256 d = np_gv.d;
00257
00258 p = np_gv.p;
00259
00260
00261 alpha=Tp->alpha;
00262
00263
00264 n=Ndoft;
00265
00266 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n\n -------------------------------------------------------------------");
00267 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n PTRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
00268 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n -------------------------------------------------------------------\n");
00269
00270
00271 Tm->updateipval ();
00272
00273
00274 capacity_matrix (lcid);
00275
00276
00277 conductivity_matrix (lcid);
00278
00279 if (istep == 0){
00280 Psolt->computation_stat (Gtt,Kmat);
00281 }
00282
00283
00284
00285 for (j=0;j<n;j++){
00286 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00287 }
00288
00289
00290
00291 Kmat->gmxv (d,p);
00292
00293
00294
00295 Kmat->scalgm (dt*alpha);
00296 Kmat->addgm (1.0,*Cmat);
00297
00298
00299 trfel_right_hand_side (0,rhs,n);
00300
00301
00302
00303 for (j=0;j<n;j++){
00304 f[j] = rhs[j] - p[j];
00305 }
00306
00307
00308
00309
00310 Psolt->par_linear_solver (Gtt,Kmat,tdlhs,f,Outt,Mesprt);
00311
00312
00313
00314 for (j=0;j<n;j++){
00315 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00316 }
00317
00318 solution_correction ();
00319
00320 approximation ();
00321
00322
00323 if (Tp->nvs==1 && Tp->pnvs==1)
00324 actual_previous_nodval ();
00325
00326 compute_req_valt (0);
00327 print_stept(0,istep,Tp->time,rhs);
00328 print_flusht();
00329
00330 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00331 solvert_save (lhs,tdlhs,f,istep,Tp->time,dt,Tp->timecont,n);
00332
00333 if (j == 0)
00334 return 0;
00335
00336 return 1;
00337 }
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 long par_one_step_nonlinear (long lcid,double time, double dt, long istep, long li, np_glob_vec &np_gv)
00349 {
00350 long j,k,n,ini;
00351 double alpha,*f,*d,*p,*lhs,*tdlhs,*rhs;
00352 double *fb,*fi;
00353 double norf_last;
00354 double zero,norf,*err,*thresh;
00355
00356
00357 Tp->time = time;
00358
00359 Tp->istep = istep;
00360
00361
00362 n=Ndoft;
00363
00364
00365 lhs = np_gv.lhs;
00366
00367 tdlhs = np_gv.tdlhs;
00368
00369 rhs = np_gv.rhs;
00370
00371 f = np_gv.f;
00372
00373 d = np_gv.d;
00374
00375 p = np_gv.p;
00376
00377
00378 fb = np_gv.fb;
00379 fi = np_gv.fi;
00380
00381
00382 nullv (fb,n);
00383 nullv (fi,n);
00384
00385
00386
00387 ini = Tp->nii;
00388
00389 err = Tp->errarr;
00390
00391 thresh=Tp->threshrhs;
00392
00393
00394 alpha=Tp->alpha;
00395
00396 zero=Tp->zero;
00397
00398 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n\n -------------------------------------------------------------------");
00399 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n PTRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
00400 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n -------------------------------------------------------------------\n");
00401
00402
00403 n=Ndoft;
00404
00405
00406 Tm->updateipval ();
00407
00408
00409 capacity_matrix (lcid);
00410
00411
00412 conductivity_matrix (lcid);
00413
00414 if (istep == 0){
00415 Psolt->computation_stat (Gtt,Kmat);
00416 }
00417
00418
00419
00420 for (j=0;j<n;j++){
00421 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00422 }
00423
00424
00425
00426 Kmat->gmxv (d,p);
00427
00428
00429
00430 Kmat->scalgm (dt*alpha);
00431 Kmat->addgm (1.0,*Cmat);
00432 Kmat->copygm (*Cmat);
00433
00434
00435 trfel_right_hand_side (0,rhs,n);
00436
00437
00438
00439 for (j=0;j<n;j++){
00440 f[j] = rhs[j] - p[j];
00441 }
00442
00443
00444
00445
00446 Psolt->par_linear_solver (Gtt,Kmat,tdlhs,f,Outt,Mesprt);
00447
00448
00449
00450 for (j=0;j<n;j++){
00451 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00452 }
00453
00454
00455 solution_correction ();
00456
00457 approximation ();
00458
00459
00460 if (Tp->nvs==1 && Tp->pnvs==1)
00461 actual_previous_nodval ();
00462
00463
00464 norf_last = 1.0e20;
00465
00466 for (j=0;j<ini;j++){
00467
00468
00469 if (Tp->trsolv == fullnewtont){
00470
00471
00472 capacity_matrix (0);
00473
00474
00475 conductivity_matrix (0);
00476
00477
00478 Kmat->gmxv (d,p);
00479
00480
00481
00482 Kmat->scalgm (dt*alpha);
00483 Kmat->addgm (1.0,*Cmat);
00484 Kmat->copygm (*Cmat);
00485 }
00486
00487 if (Tp->trestype==lrhst){
00488
00489 trfel_right_hand_side (lcid,rhs,n);
00490
00491 Kmat->gmxv (tdlhs,fi);
00492
00493 for (k=0;k<n;k++){
00494 fb[k] = rhs[k] - p[k] - fi[k];
00495 }
00496 }
00497
00498 if (Tp->trestype==fluxest){
00499
00500 internal_fluxes (fi,n);
00501
00502 for (k=0;k<n;k++){
00503 fb[k]=fi[k];
00504 }
00505 }
00506
00507 norf = Psolt->unbalanced_values (fb,Outt);
00508
00509 if (Myrank==0 && Mesprt==1) fprintf (stdout,"\n inner iteration %ld error %e",j,norf);
00510 if (norf<err[0]) break;
00511
00512 Psolt->par_linear_solver (Gtt,Cmat,fi,fb,Outt,Mesprt);
00513
00514 for (k=0;k<n;k++){
00515 tdlhs[k]+=fi[k];
00516 lhs[k]+=alpha*dt*fi[k];
00517 }
00518
00519
00520 solution_correction ();
00521
00522 approximation ();
00523
00524
00525
00526 if (Tp->nvs==1 && Tp->pnvs==1)
00527 actual_previous_nodval ();
00528
00529
00530
00531 if (Tp->convergcontrolt==yes){
00532 if (norf > norf_last){
00533
00534 if (Myrank==0 && Mesprt==1)
00535 fprintf (stdout,"\n\n convergence control: inner iteration skiped %ld error %e\n\n",j,norf);
00536
00537 return -2;
00538 }
00539
00540 norf_last = norf;
00541 }
00542 }
00543
00544 if(norf >= err[0])
00545 return -1;
00546
00547
00548
00549 compute_req_valt (lcid);
00550 print_stept(lcid,istep,Tp->time,rhs);
00551 print_flusht();
00552
00553 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00554 solvert_save (lhs,tdlhs,f,istep,Tp->time,dt,Tp->timecont,n);
00555
00556 return 0;
00557 }
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568 void par_linear_nonstat_solv_vform ()
00569 {
00570 long i,j,n;
00571 double s,zero,alpha,*d,*p,*lhs,*lhsi,*tdlhs,*rhs;
00572
00573 double dt,end_time;
00574
00575
00576 n=Ndoft;
00577
00578 lhs = Lsrst->give_lhs (0);
00579 tdlhs = Lsrst->give_tdlhs (0);
00580 lhsi = Lsrst->lhsi;
00581 rhs = Lsrst->give_rhs (0);
00582
00583 d = new double [n];
00584 p = new double [n];
00585
00586
00587
00588 nullv (lhs,n);
00589 nullv (tdlhs,n);
00590 nullv (d,n);
00591 nullv (p,n);
00592
00593
00594
00595
00596 alpha=Tp->alpha;
00597
00598 zero=Tp->zero;
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608 i=0;
00609
00610
00611 Tp->time = Tp->timecont.starttime ();
00612
00613 dt = Tp->timecont.initialtimeincr ();
00614
00615 end_time = Tp->timecont.endtime ();
00616
00617
00618
00619 approximation ();
00620 actual_previous_change ();
00621
00622 Tm->initmaterialmodels();
00623
00624
00625 print_initt(-1, "wt",Ptp->fni,Ptp->fei);
00626 print_stept(0,i,Tp->time,NULL);
00627 print_flusht();
00628
00629
00630 do{
00631
00632
00633 if (Myrank==0) fprintf (stdout,"\n iteration number %ld",i);
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 conductivity_matrix (0);
00648
00649 capacity_matrix (0);
00650
00651 if (i==0){
00652 Psolt->computation_stat (Gtt,Kmat);
00653 }
00654
00655
00656
00657
00658
00659
00660
00661 for (j=0;j<n;j++){
00662 p[j]=lhs[j]+(1.0-alpha)*dt*tdlhs[j];
00663 }
00664
00665
00666 Kmat->gmxv (p,d);
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677 Kmat->scalgm (alpha*dt);
00678 Kmat->addgm (1.0,*Cmat);
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697 trfel_right_hand_side (0,rhs,n);
00698
00699
00700 for (long ijk=0;ijk<n;ijk++){
00701 fprintf (Outt,"\n rhs %ld %le",ijk,rhs[ijk]);
00702 }
00703
00704
00705
00706 for (j=0;j<n;j++){
00707 rhs[j] = rhs[j] - d[j];
00708 d[j]=tdlhs[j];
00709 }
00710
00711
00712
00713
00714 for (long ijk=0;ijk<n;ijk++){
00715 fprintf (Outt,"\n rhs2 %ld %le",ijk,rhs[ijk]);
00716 }
00717
00718
00719
00720 Kmat->printmat (Outt);
00721
00722
00723
00724
00725
00726 Psolt->par_linear_solver (Gtt,Kmat,tdlhs,rhs,Outt,Mesprt);
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741 for (j=0;j<n;j++){
00742 s=(1.0-alpha)*d[j]+alpha*tdlhs[j];
00743 lhs[j]+=dt*s;
00744 }
00745
00746
00747
00748 for (long ijk=0;ijk<n;ijk++){
00749 fprintf (Outt,"\n lhs %ld %le",ijk,lhs[ijk]);
00750 }
00751
00752
00753 approximation ();
00754
00755
00756
00757 Tp->time = Tp->timecont.newtime ();
00758 dt = Tp->timecont.actualbacktimeincr ();
00759 i++;
00760
00761 print_stept(0,i,Tp->time,NULL);
00762 print_flusht();
00763
00764 printf ("\n Tp->time %e",Tp->time);
00765 printf ("\n i %ld",i);
00766
00767 }while(Tp->time<end_time);
00768
00769
00770
00771 delete [] p; delete [] d;
00772
00773
00774 print_closet();
00775
00776 Psolt->computation_stat_print (Outt);
00777
00778
00779 }
00780