00001 #include "backupsolt.h"
00002 #include "pnnpsolvert.h"
00003 #include "pglobalt.h"
00004 #include "seqfilest.h"
00005 #include <string.h>
00006 #include <math.h>
00007 #include "mpi.h"
00008
00009
00010
00011
00012 void par_solve_nonlinear_nonstationary_problem_dform ()
00013 {
00014 long i,j,k,n,ani,ini,lcid,stop,nsts;
00015 double dt,end_time,alpha,zero,dtmin,dtmax;
00016 double *d,*fb,*fi,*p,*f,*v,*z,*lhs,*lhsi,*lhsb,*tdlhs,*tdlhsb,*rhs,*arhs,*gz,*grhs,*err,*thresh;
00017
00018
00019 lcid=0;
00020 zero=1.0e-20;
00021
00022
00023 n=Ndoft;
00024
00025
00026 lhs = Lsrst->give_lhs (0);
00027
00028 tdlhs = Lsrst->give_tdlhs (0);
00029
00030 lhsi = Lsrst->lhsi;
00031
00032 rhs = Lsrst->give_rhs (0);
00033
00034
00035 f = new double [n];
00036 nullv (f,n);
00037
00038 d = new double [n];
00039 nullv (d,n);
00040
00041 p = new double [n];
00042 nullv (p,n);
00043
00044 v = new double [n];
00045 nullv (v,n);
00046
00047 z = new double [n];
00048 nullv (z,n);
00049
00050 lhsb = new double [n];
00051 nullv (lhsb,n);
00052
00053 tdlhsb = new double [n];
00054 nullv (tdlhs,n);
00055
00056 gz = new double [Psolt->schcom->ndofcp];
00057 grhs = new double [Psolt->schcom->ndofcp];
00058
00059
00060
00061 approximation ();
00062
00063 if (Tp->nvs==1 && Tp->pnvs==1){
00064
00065
00066 actual_previous_nodval ();
00067 }
00068
00069
00070 alpha=Tp->alpha;
00071
00072 err=Tp->errarr;
00073
00074 ini=Tp->nii;
00075
00076 thresh=Tp->threshrhs;
00077
00078
00079 dtmin=Tp->timecont.dtmin;
00080
00081 dtmax=Tp->timecont.dtmax;
00082
00083
00084 Tp->time=Tp->timecont.starttime ();
00085
00086 end_time = Tp->timecont.endtime ();
00087
00088 dt = Tp->timecont.initialtimeincr ();
00089
00090
00091
00092 Tm->initmaterialmodels();
00093
00094 approximation ();
00095
00096
00097
00098 if (Tp->nvs==1 && Tp->pnvs==1)
00099 actual_previous_nodval ();
00100
00101
00102
00103 double *uf,*buff;
00104 buff = new double [10];
00105 uf = new double [10];
00106
00107
00108
00109
00110
00111
00112
00113 i=0;
00114
00115 nsts=0;
00116
00117 stop=0;
00118
00119
00120
00121 if (Tp->hdbcont.restore_stat()){
00122 solvert_restore (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
00123 print_initt(-1, "at",Ptp->fni,Ptp->fei);
00124 print_stept(0,i,Tp->time,NULL);
00125 }
00126 else{
00127 compute_req_valt (0);
00128 print_initt(-1, "wt",Ptp->fni,Ptp->fei);
00129 print_stept(0,i,Tp->time,rhs);
00130 }
00131 print_flusht();
00132
00133
00134 Tm->initmaterialmodels();
00135
00136 do{
00137
00138
00139 Tp->time=Tp->timecont.newtime (dt);
00140
00141 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00142
00143 Tm->updateipval ();
00144 i++;
00145
00146
00147 for (j=0;j<n;j++){
00148 lhsb[j]=lhs[j];
00149 tdlhsb[j]=tdlhs[j];
00150 }
00151
00152
00153 capacity_matrix (0);
00154
00155
00156 conductivity_matrix (0);
00157
00158
00159 for (j=0;j<n;j++){
00160 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00161 }
00162
00163 fprintf (Outt,"\n\n\n\n kontrola vektoru lhs a tdlhs a d\n");
00164 for (long ijk=0;ijk<Gtt->nn;ijk++){
00165 long cn1,cn2;
00166 fprintf (Outt,"\n %10.6lf %10.6lf ",Gtt->gnodes[ijk].x,Gtt->gnodes[ijk].y);
00167 cn1=Gtt->gnodes[ijk].cn[0];
00168 cn2=Gtt->gnodes[ijk].cn[1];
00169 if (cn1>0 && cn2>0){
00170 fprintf (Outt," lhs % 12.9le tdlhs % 12.9le d % 12.9le",lhs[cn1-1],tdlhs[cn1-1],d[cn1-1]);
00171 fprintf (Outt," lhs % 12.9le tdlhs % 12.9le d % 12.9le",lhs[cn2-1],tdlhs[cn2-1],d[cn2-1]);
00172 }
00173 }
00174
00175
00176
00177 Cmat->gmxv (d,p);
00178
00179
00180
00181 Kmat->scalgm (dt*alpha);
00182 Kmat->addgm (1.0,*Cmat);
00183
00184
00185 trfel_right_hand_side (lcid,f,n);
00186
00187
00188 for (j=0;j<n;j++){
00189 rhs[j] = f[j]*alpha*dt + p[j];
00190 }
00191
00192
00193 fprintf (Outt,"\n\n\n\n kontrola vektoru p a f a rhs\n");
00194 for (long ijk=0;ijk<Gtt->nn;ijk++){
00195 long cn1,cn2;
00196 fprintf (Outt,"\n %10.6lf %10.6lf ",Gtt->gnodes[ijk].x,Gtt->gnodes[ijk].y);
00197 cn1=Gtt->gnodes[ijk].cn[0];
00198 cn2=Gtt->gnodes[ijk].cn[1];
00199 if (cn1>0 && cn2>0){
00200 fprintf (Outt," p % 12.9le f % 12.9le rhs % 12.9le",p[cn1-1],f[cn1-1],rhs[cn1-1]);
00201 fprintf (Outt," p % 12.9le f % 12.9le rhs % 12.9le",p[cn2-1],f[cn2-1],rhs[cn2-1]);
00202 }
00203 }
00204
00205
00206 nullv (z,n);
00207
00208 Kmat->diag_check (zero,rhs);
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 fprintf (Outt,"\n\n\n\n druha kontrola vektoru rhs\n");
00259 for (long ijk=0;ijk<Gtt->nn;ijk++){
00260 long cn1,cn2;
00261 fprintf (Outt,"\n %10.6lf %10.6lf ",Gtt->gnodes[ijk].x,Gtt->gnodes[ijk].y);
00262 cn1=Gtt->gnodes[ijk].cn[0];
00263 cn2=Gtt->gnodes[ijk].cn[1];
00264 if (cn1>0 && cn2>0){
00265 fprintf (Outt," rhs % 12.9le % 12.9le",rhs[cn1-1],rhs[cn2-1]);
00266 }
00267 }
00268
00269
00270
00271 Psolt->par_linear_solver (Gtt,Kmat,lhs,rhs,Outt,Mesprt);
00272
00273
00274 fprintf (Outt,"\n\n\n\n\n kontrola reseni");
00275
00276 for (long ijk=0;ijk<Gtt->nn;ijk++){
00277 long cn1,cn2;
00278 cn1=Gtt->gnodes[ijk].cn[0];
00279 cn2=Gtt->gnodes[ijk].cn[1];
00280
00281 fprintf (Outt,"\n %10.6lf %10.6lf ",Gtt->gnodes[ijk].x,Gtt->gnodes[ijk].y);
00282
00283
00284
00285
00286 if (cn1>0){
00287 fprintf (Outt," % 12.9le",lhs[cn1-1]);
00288 }
00289 else{
00290 fprintf (Outt," % 12.9le",0.0);
00291 }
00292 if (cn2>0){
00293 fprintf (Outt," % 12.9le",lhs[cn2-1]);
00294 }
00295 else{
00296 fprintf (Outt," % 12.9le",0.0);
00297 }
00298 }
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 for (j=0;j<n;j++){
00310 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
00311 }
00312
00313 solution_correction ();
00314
00315 approximation ();
00316
00317 if (Tp->nvs==1 && Tp->pnvs==1)
00318 actual_previous_nodval ();
00319
00320
00321
00322 for (j=0;j<ini;j++){
00323
00324
00325 capacity_matrix (lcid);
00326
00327
00328 conductivity_matrix (lcid);
00329
00330
00331
00332
00333 Kmat->scalgm (dt*alpha);
00334 Kmat->addgm (1.0,*Cmat);
00335
00336
00337 nullv (v,n);
00338 Kmat->gmxv (lhs,v);
00339
00340
00341 nullv (p,n);
00342 Cmat->gmxv (d,p);
00343
00344
00345
00346
00347
00348
00349 fprintf (Outt,"\n\n kontrola vektoru \n");
00350 for (k=0;k<n;k++){
00351 rhs[k]=f[k]*alpha*dt+p[k]-v[k];
00352
00353
00354
00355
00356
00357 z[k]+=f[k]*alpha*dt+p[k];
00358
00359 fprintf (Outt,"\n %6ld f,p,v,z,rhs %12.9le %12.9le %12.9le %12.9le %12.9le",k,f[k],p[k],v[k],z[k],rhs[k]);
00360
00361
00362
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398 fprintf (Outt,"\n\n kontrola vektoru rhs a z\n");
00399 for (long ijk=0;ijk<Gtt->nn;ijk++){
00400 long cn1,cn2;
00401 cn1=Gtt->gnodes[ijk].cn[0];
00402 cn2=Gtt->gnodes[ijk].cn[1];
00403 fprintf (Outt,"\n %10.6lf %10.6lf ",Gtt->gnodes[ijk].x,Gtt->gnodes[ijk].y);
00404 if (cn1>0 && cn2>0){
00405 fprintf (Outt," rhs % 12.9le z % 12.9le rhs % 12.9le z % 12.9le",rhs[cn1-1],z[cn1-1],rhs[cn2-1],z[cn2-1]);
00406 }
00407 }
00408
00409
00410
00411
00412 stop=0;
00413 for (k=0;k<Tp->ntm;k++){
00414 if (Myrank==0){
00415
00416 nullv (grhs,Psolt->schcom->ndofcp);
00417 nullv (gz,Psolt->schcom->ndofcp);
00418 }
00419
00420 stop+= Psolt->selected_norm_calculation (k,err[k],thresh[k],Gtt,rhs,grhs,z,gz,n,Outt);
00421 fprintf (stdout,"\n i,j,k %ld %ld %ld %ld",i,j,k,stop);
00422 }
00423
00424 if (stop==Tp->ntm){
00425 break;
00426 }
00427
00428 if (Mesprt != 0) fprintf (stdout,"\n iteration number %ld",j);
00429
00430 Kmat->diag_check (zero,rhs);
00431
00432
00433 Psolt->par_linear_solver (Gtt,Kmat,z,rhs,Outt,Mesprt);
00434
00435 for (k=0;k<n;k++){
00436 lhs[k]+=z[k];
00437 }
00438
00439 nullv (z,n);
00440
00441
00442 for (k=0;k<n;k++){
00443 tdlhs[k]=(lhs[k]-d[k])/dt/alpha;
00444 }
00445
00446
00447 solution_correction ();
00448
00449 approximation ();
00450
00451
00452
00453 if (Tp->nvs==1 && Tp->pnvs==1)
00454 actual_previous_nodval ();
00455
00456
00457 }
00458
00459
00460
00461 ani=j;
00462
00463 if (ani==0)
00464 nsts++;
00465
00466 if (ani==ini){
00467
00468 for (j=0;j<n;j++){
00469 lhs[j]=lhsb[j];
00470 tdlhs[j]=tdlhsb[j];
00471 }
00472
00473
00474
00475
00476
00477 solution_correction ();
00478
00479 approximation ();
00480
00481
00482
00483
00484 if (Tp->nvs==1 && Tp->pnvs==1)
00485 actual_previous_nodval ();
00486
00487
00488
00489
00490
00491 dt/=2.0;
00492 i--;
00493 Tp->timecont.oldtime ();
00494
00495
00496
00497
00498 if (dt<dtmin){
00499 fprintf (stderr,"\n\n time increment is less than minimum time increment");
00500 fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
00501 abort ();
00502 }
00503 }else{
00504
00505
00506
00507
00508 if (nsts==2){
00509 double times, acttimes;
00510 times = Tp->timecont.starttime();
00511 acttimes = Tp->timecont.actualtime();
00512
00513 if ((times+18000)>acttimes) {
00514 dt = dt;
00515 }
00516 else{
00517 dt*=2;
00518 if (dt>dtmax) {
00519 dt = dtmax;
00520 }
00521 }
00522
00523 nsts=0;
00524
00525
00526
00527
00528 }
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541 print_stept(0,i,Tp->time,NULL);
00542 print_flusht();
00543
00544 }
00545
00546 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00547 solvert_save (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
00548
00549 }while(Tp->time<end_time);
00550
00551
00552
00553
00554
00555
00556
00557 }
00558
00559
00560
00561 void par_solve_nonlinear_nonstationary_problem ()
00562 {
00563 long i,j,k,n,ini;
00564 double dt,end_time,alpha,norfb,err,totnorfb;
00565 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*arhs;
00566 MPI_Status stat;
00567
00568 n=Ndoft;
00569
00570 lhs = Lsrst->give_lhs (0);
00571 tdlhs = Lsrst->give_tdlhs (0);
00572 lhsi = Lsrst->lhsi;
00573 rhs = Lsrst->give_rhs (0);
00574
00575 d = new double [n];
00576 arhs = new double [n];
00577 p = new double [n];
00578 fb = new double [n];
00579 fi = new double [n];
00580
00581
00582 nullv (lhs,n);
00583 nullv (tdlhs,n);
00584 nullv (rhs,n);
00585 nullv (arhs,n);
00586 nullv (d,n);
00587 nullv (p,n);
00588 nullv (fb,n);
00589 nullv (fi,n);
00590
00591
00592 approximation ();
00593
00594 alpha=Tp->alpha;
00595 err=Tp->err;
00596 ini=Tp->nii;
00597
00598
00599 Tp->time=Tp->timecont.starttime ();
00600
00601 end_time = Tp->timecont.endtime ();
00602
00603 dt = Tp->timecont.initialtimeincr ();
00604
00605 print_initt(-1, "wt",Ptp->fni,Ptp->fei);
00606 print_stept(0,i,Tp->time,NULL);
00607 print_flusht();
00608
00609 long nidof=Psolt->schcom[0].nidof;
00610 long nrdof=Psolt->schcom[0].nbdof;
00611 double *uf,*buff;
00612 buff = new double [10];
00613 uf = new double [10];
00614
00615
00616
00617
00618 i=0;
00619 do{
00620
00621 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631 capacity_matrix (0);
00632
00633
00634 conductivity_matrix (0);
00635
00636 for (j=0;j<n;j++){
00637 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00638 }
00639
00640
00641 trfel_right_hand_side (0,rhs,n);
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651 Kmat->gmxv (d,p);
00652
00653
00654
00655 Kmat->scalgm (dt*alpha);
00656 Kmat->addgm (1.0,*Cmat);
00657
00658
00659 for (j=0;j<n;j++){
00660 fb[j] = rhs[j] - p[j];
00661 }
00662
00663
00664 Psolt->par_linear_solver (Gtt,Kmat,tdlhs,fb,Outt,Mesprt);
00665
00666
00667 for (j=0;j<n;j++){
00668 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00669 }
00670
00671
00672
00673
00674
00675
00676 approximation ();
00677 internal_fluxes (fi,n);
00678
00679
00680 for (j=0;j<n;j++){
00681 fb[j]=fi[j];
00682
00683 }
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757 totnorfb = Psolt->unbalanced_values (fb,Outt);
00758
00759
00760
00761
00762 if (Mesprt==1) fprintf (stdout,"\n%e %e",totnorfb,err);
00763
00764
00765 if (totnorfb<err){
00766
00767 Tp->time = Tp->timecont.newtime ();
00768 dt = Tp->timecont.actualbacktimeincr ();
00769 i++;
00770
00771 print_stept(0,i,Tp->time,NULL);
00772 print_flusht();
00773 continue;
00774 }
00775
00776
00777 for (j=0;j<ini;j++){
00778
00779
00780
00781
00782
00783
00784
00785
00786 Psolt->par_linear_solver (Gtt,Kmat,p,fb,Outt,Mesprt);
00787
00788 for (k=0;k<n;k++){
00789 tdlhs[k]+=p[k];
00790 lhs[k]+=alpha*dt*p[k];
00791 }
00792
00793
00794
00795
00796
00797
00798 approximation ();
00799 internal_fluxes (fi,n);
00800
00801
00802 for (k=0;k<n;k++){
00803 fb[k]= fi[k];
00804 }
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871 totnorfb = Psolt->unbalanced_values (fb,Outt);
00872
00873 if (Mesprt==1) fprintf (stdout,"\n iteration %ld error %e",j,totnorfb);
00874
00875 if (totnorfb<err){
00876 break;
00877 }
00878 }
00879
00880
00881
00882 Tp->time = Tp->timecont.newtime ();
00883 dt = Tp->timecont.actualbacktimeincr ();
00884 i++;
00885
00886 print_stept(0,i,Tp->time,NULL);
00887 print_flusht();
00888
00889 printf ("\n Tp->time %e",Tp->time);
00890 printf ("\n i %ld",i);
00891
00892 }while(Tp->time<end_time);
00893
00894 delete [] fi; delete [] fb; delete [] p;
00895 delete [] arhs; delete [] d;
00896
00897
00898 }
00899