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