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