00001 #include "npglvec.h"
00002 #include "npsolvert.h"
00003 #include "globalt.h"
00004 #include "globmatt.h"
00005 #include "elemswitcht.h"
00006 #include "transprint.h"
00007 #include "saddpoint.h"
00008 #include "backupsolt.h"
00009 #include <string.h>
00010 #include <errno.h>
00011
00012
00013
00014
00015
00016
00017
00018
00019 void solve_nonstationary_problem ()
00020 {
00021 if (Tp->tprob == discont_nonstat_problem){
00022
00023 linear_nonstat_solv_dform ();
00024 }
00025 else{
00026
00027
00028 nonstat_solv_vform_comp ();
00029
00030
00031
00032 }
00033 }
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 void 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
00056 approximation();
00057 nonstat_solver_init(lcid, np_gv);
00058
00059
00060 dt = Tp->timecont.initialtimeincr ();
00061
00062 dtmin=Tp->timecont.dtmin;
00063
00064 dtmax=Tp->timecont.dtmax;
00065
00066 end_time = Tp->timecont.endtime ();
00067
00068
00069 li = i = np_gv.istep;
00070
00071 nsts=0;
00072
00073 ret = 0;
00074
00075
00076 if(Tp->timecont.tct > 0){
00077 n=Ndoft;
00078 lhsb = new double [n];
00079 tdlhsb = new double [n];
00080 nullv (lhsb,n);
00081 nullv (tdlhsb,n);
00082 }
00083
00084 bool breakloop = false;
00085
00086
00087
00088 do{
00089 i++;
00090
00091 newtime = Tp->time = Tp->timecont.newtime(dt);
00092
00093 if((Tp->timecont.timefun.tfunc != stat) && Tp->timecont.tct == 0)
00094 dt = Tp->timecont.actualbacktimeincr ();
00095
00096 if(Tp->tprob == nonstationary_problem)
00097
00098
00099 ret = one_step_linear(lcid, newtime, dt, i, li, np_gv);
00100
00101 if(Tp->tprob == nonlinear_nonstationary_problem){
00102
00103
00104 ret = one_step_nonlinear(lcid, newtime, dt, i, li, np_gv);
00105 }
00106
00107 if (ret >= 0)
00108 {
00109
00110
00111 compute_req_valt (lcid);
00112 print_stept(lcid,Tp->istep,Tp->time,np_gv.rhs);
00113 print_flusht();
00114
00115 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00116 {
00117 if (Mesprt==1)
00118 fprintf (stdout,"\n Creating TRFEL backup file\n");
00119 solvert_save (np_gv.lhs,np_gv.tdlhs,np_gv.f,Tp->istep,Tp->time,dt,Tp->timecont,Ndoft);
00120 }
00121 }
00122
00123
00124 if(Tp->tprob == nonlinear_nonstationary_problem && Tp->timecont.tct > 0)
00125 {
00126
00127 if (ret >= 0)
00128 {
00129 for (k=0;k<n;k++){
00130 lhsb[k] = np_gv.lhs[k];
00131 tdlhsb[k] = np_gv.tdlhs[k];
00132 }
00133
00134 approximation ();
00135
00136 if (ret == 0)
00137 nsts++;
00138 dtdef = Tp->timecont.actualforwtimeincr();
00139 if (nsts==2)
00140 {
00141 dt*=2.0;
00142 nsts=0;
00143
00144 if (Mesprt==1)
00145 fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00146 }
00147 if (dt<=dtdef)
00148 dt = dtdef;
00149
00150 if (dt>dtmax)
00151 dt = dtmax;
00152 }
00153 else
00154 {
00155
00156
00157 dt/=2.0;
00158 Tp->timecont.oldtime ();
00159 i--;
00160
00161 for (k=0;k<n;k++){
00162 np_gv.lhs[k] = lhsb[k];
00163 np_gv.tdlhs[k] = tdlhsb[k];
00164 }
00165
00166
00167 approximation ();
00168
00169 if (Mesprt==1)
00170 fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00171
00172 if (dt<dtmin)
00173 {
00174 if (Mesprt==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00175 if (Mesprt==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00176 if (Mesprt==1) fprintf (stderr,"\n FORCED output of results from this step is performed\n");
00177 compute_req_valt (lcid);
00178 print_stept_forced(lcid, i, Tp->time, np_gv.rhs);
00179 print_flusht();
00180 break;
00181 }
00182 }
00183 }
00184
00185
00186 if (Tp->adaptivityflag && Tp->time < end_time)
00187 if (i%2 == 0)
00188 breakloop = Adat->run (2, true);
00189
00190 }while(Tp->time < end_time && breakloop == false);
00191
00192 print_closet();
00193
00194 if(Tp->timecont.tct > 0){
00195 delete [] lhsb;
00196 delete [] tdlhsb;
00197 }
00198 }
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208 void nonstat_solver_init (long lcid, np_glob_vec &np_gv)
00209 {
00210 long n;
00211 double dt;
00212
00213
00214 n=Ndoft;
00215
00216 np_gv.alloc(n);
00217
00218 np_gv.lhs = Lsrst->give_lhs (lcid);
00219
00220 np_gv.tdlhs = Lsrst->give_tdlhs (lcid);
00221
00222 np_gv.rhs = Lsrst->give_rhs (lcid);
00223
00224 if(Tp->tprob == nonlinear_nonstationary_problem)
00225 np_gv.alloc_aux(n);
00226
00227
00228 Tp->time=Tp->timecont.starttime ();
00229
00230 dt = Tp->timecont.initialtimeincr ();
00231
00232
00233 if (Tp->nvs==1 && Tp->pnvs==1)
00234 actual_previous_nodval ();
00235
00236
00237 if ((Adat != NULL) && (Adat->give_step()>0)) {
00238 Adat->statedata_restore();
00239
00240
00241 print_initt(-1, "at");
00242 }
00243 else {
00244
00245 if (Tp->hdbcont.restore_stat()){
00246 if (Mesprt==1)
00247 fprintf (stdout,"\n Reading of TRFEL backup file\n");
00248 solvert_restore (np_gv.lhs, np_gv.tdlhs, np_gv.f, np_gv.istep, Tp->time, dt, Tp->timecont, n);
00249
00250 Tm->initmaterialmodels();
00251 compute_req_valt (lcid);
00252 print_initt(-1, "at");
00253 }
00254 else{
00255
00256 Tm->initmaterialmodels();
00257 compute_req_valt (lcid);
00258 print_initt(-1, "wt");
00259 print_stept(lcid, np_gv.istep, Tp->time, np_gv.rhs);
00260 }
00261 print_flusht();
00262 }
00263
00264
00265 if (Adat != NULL){
00266 if (Adat->give_step()) {
00267 Adat->run (2, true);
00268 Adat->answer = 0;
00269 }
00270 }
00271 }
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284 long one_step_linear (long lcid,double time, double dt, long istep, long li, np_glob_vec &np_gv)
00285 {
00286 long j,n;
00287 double alpha,*f,*d,*p,*lhs,*tdlhs,*rhs;
00288
00289
00290 Tp->time = time;
00291
00292 Tp->istep = istep;
00293 Tp->jstep = -1;
00294
00295
00296 lhs = np_gv.lhs;
00297
00298 tdlhs = np_gv.tdlhs;
00299
00300 rhs = np_gv.rhs;
00301
00302 f = np_gv.f;
00303
00304 d = np_gv.d;
00305
00306 p = np_gv.p;
00307
00308
00309 alpha=Tp->alpha;
00310
00311
00312 n=Ndoft;
00313
00314 if (Mesprt==1){
00315 fprintf (stdout,"\n\n ------------------------------------------------------------------------");
00316 fprintf (stdout,"\n TRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
00317 fprintf (stdout,"\n ------------------------------------------------------------------------\n");
00318 }
00319
00320
00321 Tm->updateipval ();
00322
00323
00324 capacity_matrix (lcid);
00325
00326
00327 conductivity_matrix (lcid);
00328
00329
00330
00331 for (j=0;j<n;j++){
00332 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00333 }
00334
00335
00336
00337 Kmat->gmxv (d,p);
00338
00339
00340
00341 Kmat->scalgm (dt*alpha);
00342 Kmat->addgm (1.0,*Cmat);
00343
00344
00345 trfel_right_hand_side (0,rhs,n);
00346
00347
00348
00349 for (j=0;j<n;j++){
00350 f[j] = rhs[j] - p[j];
00351 }
00352
00353
00354
00355
00356 Tp->ssle->solve_system (Gtt,Kmat,tdlhs,f,Outt);
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379 for (j=0;j<n;j++){
00380 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00381 }
00382
00383 solution_correction ();
00384
00385 approximation ();
00386
00387
00388 if (Tp->nvs==1 && Tp->pnvs==1)
00389 actual_previous_nodval ();
00390
00391 return 0;
00392 }
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405 long one_step_nonlinear (long lcid,double time, double dt, long istep, long li, np_glob_vec &np_gv)
00406 {
00407 long j,k,n,ini;
00408 double alpha,*f,*d,*p,*lhs,*tdlhs,*rhs;
00409 double *fb,*fi;
00410 double norf_last;
00411 double zero,norf,*err,*thresh;
00412
00413
00414 Tp->time = time;
00415
00416 Tp->istep = istep;
00417 Tp->jstep = -1;
00418
00419
00420 n=Ndoft;
00421
00422
00423 lhs = np_gv.lhs;
00424
00425 tdlhs = np_gv.tdlhs;
00426
00427 rhs = np_gv.rhs;
00428
00429 f = np_gv.f;
00430
00431 d = np_gv.d;
00432
00433 p = np_gv.p;
00434
00435
00436 fb = np_gv.fb;
00437 fi = np_gv.fi;
00438
00439
00440 nullv (fb,n);
00441 nullv (fi,n);
00442
00443
00444
00445 ini = Tp->nii;
00446
00447 err = Tp->errarr;
00448
00449 thresh=Tp->threshrhs;
00450
00451
00452 alpha=Tp->alpha;
00453
00454 zero=Tp->zero;
00455
00456 if (Mesprt==1) fprintf (stdout,"\n\n --------------------------------------------------------------");
00457 if (Mesprt==1) fprintf (stdout,"\n TRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
00458 if (Mesprt==1) fprintf (stdout,"\n --------------------------------------------------------------\n");
00459
00460
00461 n=Ndoft;
00462
00463
00464 Tm->updateipval ();
00465
00466
00467 capacity_matrix (lcid);
00468
00469
00470 conductivity_matrix (lcid);
00471
00472
00473
00474 for (j=0;j<n;j++){
00475 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00476 }
00477
00478
00479
00480 Kmat->gmxv (d,p);
00481
00482
00483
00484 Kmat->scalgm (dt*alpha);
00485 Kmat->addgm (1.0,*Cmat);
00486 Kmat->copygm (*Cmat);
00487
00488
00489 trfel_right_hand_side (0,rhs,n);
00490
00491
00492
00493 for (j=0;j<n;j++){
00494 f[j] = rhs[j] - p[j];
00495 }
00496
00497
00498
00499 Tp->ssle->solve_system (Gtt,Cmat,tdlhs,f,Outt);
00500
00501
00502
00503 for (j=0;j<n;j++){
00504 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00505 }
00506
00507
00508
00509 solution_correction ();
00510
00511
00512 approximation ();
00513
00514
00515 if (Tp->nvs==1 && Tp->pnvs==1)
00516 actual_previous_nodval ();
00517
00518
00519
00520 norf_last = 1.0e20;
00521
00522 for (j=0;j<ini;j++){
00523 Tp->jstep = j;
00524
00525 if (Tp->trsolv == fullnewtont){
00526
00527
00528 capacity_matrix (0);
00529
00530
00531 conductivity_matrix (0);
00532
00533
00534 Kmat->gmxv (d,p);
00535
00536
00537
00538 Kmat->scalgm (dt*alpha);
00539 Kmat->addgm (1.0,*Cmat);
00540 Kmat->copygm (*Cmat);
00541 }
00542
00543 if (Tp->trestype==lrhst){
00544
00545 trfel_right_hand_side (lcid,rhs,n);
00546
00547
00548 Kmat->gmxv (tdlhs,fi);
00549
00550 for (k=0;k<n;k++){
00551 fb[k] = rhs[k] - p[k] - fi[k];
00552 }
00553 }
00554
00555 if (Tp->trestype==fluxest){
00556
00557 internal_fluxes (fi,n);
00558
00559 for (k=0;k<n;k++){
00560 fb[k]=fi[k];
00561 }
00562 }
00563
00564 norf = ss (fb,fb,n);
00565 if (Mesprt==1) fprintf (stdout,"\n inner iteration %ld error %e",j,norf);
00566 if (norf<err[0]) break;
00567
00568 Tp->ssle->solve_system (Gtt,Cmat,fi,fb,Outt);
00569
00570 for (k=0;k<n;k++){
00571 tdlhs[k]+=fi[k];
00572 lhs[k]+=alpha*dt*fi[k];
00573 }
00574
00575
00576 solution_correction ();
00577
00578 approximation ();
00579
00580
00581
00582 if (Tp->nvs==1 && Tp->pnvs==1)
00583 actual_previous_nodval ();
00584
00585
00586
00587 if (Tp->convergcontrolt==yes){
00588 if (norf > norf_last){
00589
00590 if (Mesprt==1) fprintf (stdout,"\n\n divergence control: inner iteration skiped %ld error %e\n\n",j,norf);
00591
00592 return -2;
00593 }
00594 norf_last = norf;
00595 }
00596
00597
00598 }
00599
00600 if(norf >= err[0])
00601 return -1;
00602
00603 return 0;
00604 }
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 void linear_nonstat_solv_vform ()
00617 {
00618 long i,j,n;
00619 double dt,end_time,alpha,*d,*p,*lhs,*tdlhs,*rhs,*f;
00620
00621
00622 n=Ndoft;
00623
00624
00625 lhs = Lsrst->give_lhs (0);
00626 nullv (lhs,n);
00627
00628 tdlhs = Lsrst->give_tdlhs (0);
00629 nullv (tdlhs,n);
00630
00631 rhs = Lsrst->give_rhs (0);
00632 nullv (rhs,n);
00633
00634
00635 f = new double [n];
00636 nullv (f,n);
00637
00638 d = new double [n];
00639 nullv (d,n);
00640
00641 p = new double [n];
00642 nullv (p,n);
00643
00644
00645 alpha=Tp->alpha;
00646
00647
00648 Tp->time=Tp->timecont.starttime ();
00649
00650 end_time = Tp->timecont.endtime ();
00651
00652 dt = Tp->timecont.initialtimeincr ();
00653
00654 i=0;
00655
00656
00657 Tm->initmaterialmodels();
00658
00659 approximation();
00660
00661
00662 if (Tp->nvs==1 && Tp->pnvs==1)
00663 actual_previous_nodval ();
00664
00665
00666 if ((Adat != NULL) && (Adat->give_step()>0)) {
00667 Adat->statedata_restore();
00668
00669 i = Tp->istep;
00670 print_initt(-1, "at");
00671 }
00672 else {
00673
00674 if (Tp->hdbcont.restore_stat()){
00675 solvert_restore (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
00676 print_initt(-1, "at");
00677
00678 }
00679 else{
00680 compute_req_valt (0);
00681 print_initt(-1, "wt");
00682 print_stept(0,i,Tp->time,rhs);
00683 }
00684 }
00685
00686
00687 if (Adat != NULL){
00688 if (Adat->give_step()) {
00689 Adat->run (2, true);
00690 Adat->answer = 0;
00691 }
00692 }
00693
00694 bool breakloop = false;
00695
00696
00697
00698 do{
00699
00700
00701 Tp->time=Tp->timecont.newtime ();
00702
00703 dt = Tp->timecont.actualbacktimeincr ();
00704
00705 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00706
00707 Tm->updateipval ();
00708
00709 i++;
00710 Tp->istep = i;
00711
00712
00713 capacity_matrix (0);
00714
00715
00716
00717
00718 conductivity_matrix (0);
00719
00720
00721
00722
00723 for (j=0;j<n;j++){
00724 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00725 }
00726
00727
00728
00729 Kmat->gmxv (d,p);
00730
00731
00732
00733
00734 Kmat->scalgm (dt*alpha);
00735 Kmat->addgm (1.0,*Cmat);
00736
00737
00738 trfel_right_hand_side (0,rhs,n);
00739
00740
00741
00742 for (j=0;j<n;j++){
00743 f[j] = rhs[j] - p[j];
00744 }
00745
00746
00747
00748
00749 Tp->ssle->solve_system (Gtt,Kmat,tdlhs,f,Outt);
00750
00751
00752
00753 for (j=0;j<n;j++){
00754
00755
00756 if(lhs[j]>=0.0)
00757 lhs[j]=-1.0;
00758
00759 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00760
00761
00762
00763
00764 }
00765
00766 solution_correction ();
00767
00768 approximation ();
00769 compute_req_valt (0);
00770
00771
00772 if (Tp->nvs==1 && Tp->pnvs==1)
00773 actual_previous_nodval ();
00774
00775 print_stept(0,i,Tp->time,rhs);
00776 print_flusht();
00777
00778 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00779 {
00780 if (Mesprt==1)
00781 fprintf (stdout,"\n Creating TRFEL backup file\n");
00782 solvert_save (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
00783 }
00784
00785
00786 if (Tp->adaptivityflag && Tp->time < end_time)
00787 if (i%2 == 0)
00788 breakloop = Adat->run (2, true);
00789
00790 }while(Tp->time < end_time && breakloop == false);
00791
00792 delete [] p;
00793 delete [] d;
00794 delete [] f;
00795
00796 print_closet();
00797 }
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810 void linear_nonstat_solv_dform ()
00811 {
00812 long i,j,n;
00813 double dt,end_time,alpha;
00814 double *d,*f,*p,*lhs,*tdlhs,*rhs,*fi;
00815
00816
00817 n=Ndoft;
00818
00819
00820 lhs = Lsrst->give_lhs (0);
00821 nullv (lhs,n);
00822
00823 tdlhs = Lsrst->give_tdlhs (0);
00824 nullv (tdlhs,n);
00825
00826 rhs = Lsrst->give_rhs (0);
00827 nullv (rhs,n);
00828
00829
00830 f = new double [n];
00831 nullv (f,n);
00832
00833 d = new double [n];
00834 nullv (d,n);
00835
00836 p = new double [n];
00837 nullv (p,n);
00838
00839 fi = new double [n];
00840 nullv (fi,n);
00841
00842
00843
00844 approximation ();
00845
00846
00847
00848 if (Tp->nvs==1 && Tp->pnvs==1)
00849 actual_previous_nodval ();
00850
00851
00852
00853 alpha=Tp->alpha;
00854
00855
00856 Tp->time=Tp->timecont.starttime ();
00857
00858 end_time = Tp->timecont.endtime ();
00859
00860 dt = Tp->timecont.initialtimeincr ();
00861
00862
00863 i=0;
00864
00865
00866 Tm->initmaterialmodels();
00867
00868
00869 if (Tp->hdbcont.restore_stat()){
00870 solvert_restore (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
00871 print_initt(-1, "at");
00872
00873 }
00874 else{
00875 compute_req_valt (0);
00876 print_initt(-1, "wt");
00877 print_stept(0,i,Tp->time,rhs);
00878 }
00879
00880
00881
00882
00883 do{
00884
00885
00886 Tp->time=Tp->timecont.newtime ();
00887
00888 dt = Tp->timecont.actualbacktimeincr ();
00889
00890 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00891
00892 Tm->updateipval ();
00893
00894 i++;
00895
00896
00897 capacity_matrix (0);
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907 conductivity_matrix (0);
00908
00909
00910
00911
00912
00913
00914
00915
00916 for (j=0;j<n;j++){
00917 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00918 }
00919
00920
00921
00922 Cmat->gmxv (d,p);
00923
00924
00925
00926
00927 Kmat->scalgm (dt*alpha);
00928 Kmat->addgm (1.0,*Cmat);
00929
00930
00931
00932 trfel_right_hand_side (0,rhs,n);
00933
00934
00935
00936
00937 for (j=0;j<n;j++){
00938 f[j] = rhs[j]*alpha*dt + p[j];
00939 }
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954 Tp->ssle->solve_system (Gtt,Kmat,lhs,f,Outt);
00955
00956
00957
00958 for (j=0;j<n;j++){
00959 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
00960 }
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972 solution_correction ();
00973 approximation ();
00974
00975
00976 if (Tp->nvs==1 && Tp->pnvs==1)
00977 actual_previous_nodval ();
00978
00979
00980 if (Tp->fluxcomp >= 1)
00981 internal_fluxes (fi,n);
00982
00983 print_stept(0,i,Tp->time,rhs);
00984 print_stept(0,i,Tp->time,NULL);
00985 print_flusht();
00986
00987 if ((Tp->timecont.isitimptime ()==1) && (Tp->hdbcont.hdbtype))
00988 {
00989 if (Mesprt==1)
00990 fprintf (stdout,"\n Creating TRFEL backup file\n");
00991 solvert_save (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
00992 }
00993
00994 }while(Tp->time<end_time);
00995
00996 delete [] fi;
00997 delete [] f;
00998 delete [] p;
00999 delete [] d;
01000
01001 print_closet();
01002 }
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013 void linear_nonstat_radiation_solv_dform ()
01014 {
01015 long i,j,n;
01016 double dt,end_time,alpha;
01017 double *d,*f,*p,*lhs,*tdlhs,*rhs,*fi;
01018
01019
01020 n=Ndoft;
01021
01022
01023 lhs = Lsrst->give_lhs (0);
01024 nullv (lhs,n);
01025
01026 tdlhs = Lsrst->give_tdlhs (0);
01027 nullv (tdlhs,n);
01028
01029 rhs = Lsrst->give_rhs (0);
01030 nullv (rhs,n);
01031
01032
01033 f = new double [n];
01034 nullv (f,n);
01035
01036 d = new double [n];
01037 nullv (d,n);
01038
01039 p = new double [n];
01040 nullv (p,n);
01041
01042 fi = new double [n];
01043 nullv (fi,n);
01044
01045
01046
01047 approximation ();
01048
01049
01050
01051 if (Tp->nvs==1 && Tp->pnvs==1)
01052 actual_previous_nodval ();
01053
01054
01055
01056 alpha=Tp->alpha;
01057
01058
01059 Tp->time=Tp->timecont.starttime ();
01060
01061 end_time = Tp->timecont.endtime ();
01062
01063 dt = Tp->timecont.initialtimeincr ();
01064
01065
01066 i=0;
01067
01068
01069 Tm->initmaterialmodels();
01070
01071
01072 if (Tp->hdbcont.restore_stat()){
01073 solvert_restore (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
01074 print_initt(-1, "at");
01075
01076 }
01077 else{
01078 compute_req_valt (0);
01079 print_initt(-1, "wt");
01080 print_stept(0,i,Tp->time,rhs);
01081 }
01082
01083
01084
01085
01086 do{
01087
01088
01089 Tp->time=Tp->timecont.newtime ();
01090
01091 dt = Tp->timecont.actualbacktimeincr ();
01092
01093 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
01094
01095 Tm->updateipval ();
01096
01097 i++;
01098
01099
01100 capacity_matrix (0);
01101
01102
01103
01104 conductivity_matrix (0);
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115 for (j=0;j<n;j++){
01116 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
01117 }
01118
01119
01120
01121 Cmat->gmxv (d,p);
01122
01123
01124
01125
01126 Kmat->scalgm (dt*alpha);
01127 Kmat->addgm (1.0,*Cmat);
01128
01129
01130
01131 trfel_right_hand_side (0,rhs,n);
01132
01133
01134
01135
01136 for (j=0;j<n;j++){
01137 f[j] = rhs[j]*alpha*dt + p[j];
01138 }
01139
01140 Tt->edge_temperature ();
01141 Tt->heat_fluxes (f,Outt);
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155 Tp->ssle->solve_system (Gtt,Kmat,lhs,f,Outt);
01156
01157
01158
01159 for (j=0;j<n;j++){
01160 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
01161 }
01162
01163
01164 solution_correction ();
01165 approximation ();
01166
01167
01168 if (Tp->nvs==1 && Tp->pnvs==1)
01169 actual_previous_nodval ();
01170
01171
01172 if (Tp->fluxcomp >= 1)
01173 internal_fluxes (fi,n);
01174
01175 print_stept(0,i,Tp->time,rhs);
01176 print_stept(0,i,Tp->time,NULL);
01177 print_flusht();
01178
01179 if ((Tp->timecont.isitimptime ()==1) && (Tp->hdbcont.hdbtype))
01180 {
01181 if (Mesprt==1)
01182 fprintf (stdout,"\n Creating TRFEL backup file\n");
01183 solvert_save (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
01184 }
01185 }while(Tp->time<end_time);
01186
01187 delete [] fi;
01188 delete [] f;
01189 delete [] p;
01190 delete [] d;
01191
01192 print_closet();
01193 }
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206 void linear_nonstat_solv_dform_subcycl ()
01207 {
01208 long i,j,n;
01209 double dt,end_time,alpha;
01210 double *d,*f,*p,*lhs,*tdlhs,*rhs,*fi;
01211
01212
01213 n=Ndoft;
01214
01215
01216 lhs = Lsrst->give_lhs (0);
01217 nullv (lhs,n);
01218
01219 tdlhs = Lsrst->give_tdlhs (0);
01220 nullv (tdlhs,n);
01221
01222 rhs = Lsrst->give_rhs (0);
01223 nullv (rhs,n);
01224
01225
01226 f = new double [n];
01227 nullv (f,n);
01228
01229 d = new double [n];
01230 nullv (d,n);
01231
01232 p = new double [n];
01233 nullv (p,n);
01234
01235 fi = new double [n];
01236 nullv (fi,n);
01237
01238
01239
01240 approximation ();
01241
01242
01243
01244 if (Tp->nvs==1 && Tp->pnvs==1)
01245 actual_previous_nodval ();
01246
01247
01248
01249 alpha=Tp->alpha;
01250
01251
01252 Tp->time=Tp->timecont.starttime ();
01253
01254 end_time = Tp->timecont.endtime ();
01255
01256 dt = Tp->timecont.initialtimeincr ();
01257
01258
01259 i=0;
01260
01261
01262 Tm->initmaterialmodels();
01263
01264
01265 if (Tp->hdbcont.restore_stat()){
01266 solvert_restore (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
01267 print_initt(-1, "at");
01268
01269 }
01270 else{
01271 compute_req_valt (0);
01272 print_initt(-1, "wt");
01273 print_stept(0,i,Tp->time,rhs);
01274 }
01275
01276
01277 Gtt->stop->elem_lists ();
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288 do{
01289
01290
01291
01292 Tp->time=Tp->timecont.newtime ();
01293
01294 dt = Tp->timecont.actualbacktimeincr ();
01295
01296 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
01297
01298 Tm->updateipval ();
01299
01300 i++;
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315 capacity_matrix (0);
01316
01317
01318 conductivity_matrix (0);
01319
01320
01321
01322
01323
01324
01325
01326 capacity_matrix (0);
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336 conductivity_matrix (0);
01337
01338
01339
01340
01341
01342
01343
01344
01345 for (j=0;j<n;j++){
01346 d[j]=lhs[j] + (1.0-alpha)*dt*tdlhs[j];
01347 }
01348
01349
01350
01351 Cmat->gmxv (d,p);
01352
01353
01354
01355
01356 Kmat->scalgm (dt*alpha);
01357 Kmat->addgm (1.0,*Cmat);
01358
01359
01360
01361 trfel_right_hand_side (0,rhs,n);
01362
01363
01364
01365
01366 for (j=0;j<n;j++){
01367 f[j] = rhs[j]*alpha*dt + p[j];
01368 }
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383 Tp->ssle->solve_system (Gtt,Kmat,lhs,f,Outt);
01384
01385
01386
01387 for (j=0;j<n;j++){
01388 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
01389 }
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401 solution_correction ();
01402 approximation ();
01403
01404
01405 if (Tp->nvs==1 && Tp->pnvs==1)
01406 actual_previous_nodval ();
01407
01408
01409 if (Tp->fluxcomp >= 1)
01410 internal_fluxes (fi,n);
01411
01412 print_stept(0,i,Tp->time,rhs);
01413 print_stept(0,i,Tp->time,NULL);
01414 print_flusht();
01415
01416 if ((Tp->timecont.isitimptime ()==1) && (Tp->hdbcont.hdbtype))
01417 {
01418 if (Mesprt==1)
01419 fprintf (stdout,"\n Creating TRFEL backup file\n");
01420 solvert_save (lhs,tdlhs,f,i,Tp->time,dt,Tp->timecont,n);
01421 }
01422 }while(Tp->time<end_time);
01423
01424 delete [] fi;
01425 delete [] f;
01426 delete [] p;
01427 delete [] d;
01428
01429 print_closet();
01430 }