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