00001 #include "nspsolvert.h"
00002 #include "globalt.h"
00003 #include "globmatt.h"
00004 #include "elemswitcht.h"
00005 #include "transprint.h"
00006 #include <string.h>
00007 #include "npsolvert.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 void solve_nonlinear_stationary_problem_pokus ()
00018 {
00019 long j,k,n,lcid,ni,ini;
00020 double zero,err,lambda,dlambda,dlambdamin,dlambdamax,norf;
00021 double *r,*rb,*dr,*f,*fi,*fb,*fp,*rhs,*jk;
00022
00023
00024 lcid=0;
00025
00026 n = Ndoft;
00027
00028 ni = Tp->nlman.ninr;
00029
00030 ini = Tp->nlman.niilnr;
00031
00032 zero = Tp->zero;
00033
00034 err = Tp->nlman.errnr;
00035
00036 dlambda=Tp->nlman.incrnr;
00037
00038 dlambdamin=Tp->nlman.minincrnr;
00039
00040 dlambdamax=Tp->nlman.maxincrnr;
00041
00042 rhs=Lsrst->give_rhs (0);
00043
00044
00045
00046 rb = new double [n];
00047 f = new double [n];
00048 fb = new double [n];
00049 fi = new double [n];
00050 dr = new double [n];
00051 jk = new double [n];
00052 memset (rb,0,n*sizeof(double));
00053 memset (f,0,n*sizeof(double));
00054 memset (fb,0,n*sizeof(double));
00055 memset (fi,0,n*sizeof(double));
00056 memset (dr,0,n*sizeof(double));
00057
00058
00059 Tm->initmaterialmodels();
00060
00061 approximation ();
00062
00063
00064 trfel_right_hand_side (lcid,rhs,n);
00065
00066
00067 r = Lsrst->give_lhs (lcid);
00068
00069 fp = Lsrst->give_rhs (lcid*2);
00070
00071
00072
00073
00074
00075 lambda=0.0;
00076
00077
00078 conductivity_matrix (lcid);
00079
00080 copyv (fp,fb,n);
00081
00082
00083 Tp->ssle->solve_system (Gtt,Kmat,r,fb,Outt);
00084
00085 nullv (fb,n);
00086
00087
00088
00089
00090 print_initt(-1, "wt");
00091
00092
00093 for (j=0;j<ini;j++){
00094
00095 if(Tp->homogt != 1)
00096 solution_correction ();
00097
00098 compute_req_valt (0);
00099 approximation ();
00100
00101
00102 conductivity_matrix (0);
00103
00104 if (Tp->trestype==lrhst){
00105
00106 trfel_right_hand_side (0,rhs,n);
00107
00108 Kmat->gmxv (r,fb);
00109
00110 for (k=0;k<n;k++){
00111 fi[k] = rhs[k] - fb[k];
00112 }
00113 }
00114
00115 if (Tp->trestype==fluxest){
00116 internal_fluxes (fi,n);
00117 }
00118
00119
00120 norf=normv(fi,n);
00121
00122
00123 fprintf (stdout,"\n inner loop j %ld norf=%20.16le dlambda %e",j,norf,dlambda);
00124
00125 if (norf<err)
00126 break;
00127
00128
00129 Tp->ssle->solve_system (Gtt,Kmat,dr,fi,Outt);
00130
00131 for (k=0;k<n;k++){
00132 r[k]+=dr[k];
00133 }
00134
00135 }
00136
00137 print_stept(lcid, 0, lambda, f);
00138 print_flusht();
00139 print_closet();
00140
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 void solve_nonlinear_stationary_problem ()
00152 {
00153 long i,j,k,n,lcid,ni,ini;
00154 double zero,err,lambda,dlambda,dlambdamin,dlambdamax,norf;
00155 double *r,*rb,*dr,*f,*fi,*fb,*fp,*rhs,*lhsi,*jk;
00156
00157
00158 lcid=0;
00159
00160 n = Ndoft;
00161
00162 ni = Tp->nlman.ninr;
00163
00164 ini = Tp->nlman.niilnr;
00165
00166 zero = Tp->zero;
00167
00168 err = Tp->nlman.errnr;
00169
00170 dlambda=Tp->nlman.incrnr;
00171
00172 dlambdamin=Tp->nlman.minincrnr;
00173
00174 dlambdamax=Tp->nlman.maxincrnr;
00175 i = 0;
00176
00177 rhs=Lsrst->give_rhs (0);
00178 lhsi=Lsrst->give_lhsi (0);
00179
00180
00181 rb = new double [n];
00182 f = new double [n];
00183 fb = new double [n];
00184 fi = new double [n];
00185 dr = new double [n];
00186 jk = new double [n];
00187 memset (rb,0,n*sizeof(double));
00188 memset (f,0,n*sizeof(double));
00189 memset (fb,0,n*sizeof(double));
00190 memset (fi,0,n*sizeof(double));
00191 memset (dr,0,n*sizeof(double));
00192
00193
00194 trfel_right_hand_side (lcid,rhs,n);
00195
00196
00197 r = Lsrst->give_lhs (lcid);
00198
00199 fp = Lsrst->give_rhs (lcid*2);
00200
00201
00202
00203
00204
00205 lambda=0.0;
00206
00207
00208
00209
00210 conductivity_matrix (lcid);
00211
00212 copyv (fp,fb,n);
00213 copyv (rhs,fi,n);
00214
00215
00216 Tp->ssle->solve_system (Gtt,Kmat,r,fb,Outt);
00217
00218
00219 conductivity_matrix (lcid);
00220
00221
00222
00223
00224
00225
00226
00227
00228 print_initt(-1, "wt");
00229 print_stept(lcid, 0, 0.0, f);
00230 print_flusht();
00231
00232 solution_correction ();
00233 compute_req_valt (0);
00234 approximation ();
00235
00236
00237
00238
00239
00240
00241 norf=normv(fi,n);
00242
00243
00244 if (norf<err){
00245
00246 print_stept(lcid, i+1, lambda, f);
00247 print_flusht();
00248
00249 }
00250 else{
00251
00252
00253 for (j=0;j<ini;j++){
00254
00255 conductivity_matrix (lcid);
00256
00257 Tp->ssle->solve_system (Gtt,Kmat,dr,fi,Outt);
00258
00259 for (k=0;k<n;k++){
00260 r[k]+=dr[k];
00261 }
00262
00263 conductivity_matrix (lcid);
00264
00265
00266
00267 solution_correction ();
00268 compute_req_valt (0);
00269 approximation ();
00270
00271
00272
00273
00274 norf=normv(fi,n);
00275
00276
00277 fprintf (stdout,"\n increment %ld inner loop j %ld norf=%20.16le dlambda %e",i,j,norf,dlambda);
00278
00279
00280 if (norf<err)
00281 break;
00282 }
00283 }
00284
00285
00286 print_stept(lcid, i+1, lambda, f);
00287 print_flusht();
00288 print_closet();
00289
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299 void solve_nonlinear_stationary_problem_old ()
00300 {
00301 long i,j,k,n,lcid,ni,ini,modif;
00302 double zero,err,lambda,blambda,dlambda,dlambdamin,dlambdamax,norf,norfa;
00303 double *d,*dd,*db,*fc,*fp,*fb,*f,*rhs;
00304
00305
00306 lcid=0;
00307
00308 n = Ndoft;
00309
00310 ni = Tp->nlman.ninr;
00311
00312 ini = Tp->nlman.niilnr;
00313
00314 zero = Tp->zero;
00315
00316 err = Tp->nlman.errnr;
00317
00318 dlambda=Tp->nlman.incrnr;
00319
00320 dlambdamin=Tp->nlman.minincrnr;
00321
00322 dlambdamax=Tp->nlman.maxincrnr;
00323
00324
00325 dd = new double [n];
00326 memset (dd,0,n*sizeof(double));
00327 db = new double [n];
00328 memset (db,0,n*sizeof(double));
00329 fc = new double [n];
00330 memset (fc,0,n*sizeof(double));
00331 fp = new double [n];
00332 memset (fp,0,n*sizeof(double));
00333 fb = new double [n];
00334 memset (fb,0,n*sizeof(double));
00335 f = new double [n];
00336 memset (f,0,n*sizeof(double));
00337
00338
00339
00340
00341 rhs=Lsrst->give_rhs (0);
00342
00343
00344 trfel_right_hand_side (lcid,rhs,n);
00345
00346
00347 d = Lsrst->give_lhs (lcid);
00348
00349 fp = Lsrst->give_rhs (lcid*2);
00350
00351 fc = Lsrst->give_rhs (lcid*2+1);
00352
00353
00354 lambda=0.0;
00355
00356
00357
00358 conductivity_matrix (lcid);
00359
00360 copyv (fc,fb,n);
00361
00362
00363 Tp->ssle->solve_system (Gtt,Kmat,d,fb,Outt);
00364
00365
00366
00367
00368 print_initt(-1, "wt");
00369 print_stept(lcid, 0, 0.0, f);
00370 print_flusht();
00371
00372 for (i=0;i<ni;i++){
00373
00374
00375 for (j=0;j<n;j++){
00376 db[j]=d[j];
00377 }
00378
00379 blambda=lambda;
00380
00381 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
00382
00383
00384 for (j=0;j<n;j++){
00385 f[j]=fc[j]+(lambda+dlambda)*fp[j];
00386 fb[j]=dlambda*fp[j];
00387 }
00388
00389
00390 conductivity_matrix (lcid);
00391
00392
00393 Tp->ssle->solve_system (Gtt,Kmat,dd,fb,Outt);
00394
00395 for (j=0;j<n;j++){
00396 d[j]+=dd[j];
00397 }
00398
00399 conductivity_matrix (lcid);
00400
00401 solution_correction ();
00402 compute_req_valt (0);
00403 approximation ();
00404
00405
00406
00407 Kmat->gmxv (d,fb);
00408
00409
00410 for (j=0;j<n;j++){
00411 fb[j]=f[j]-fb[j];
00412 }
00413
00414 norf=normv(fb,n);
00415
00416 if (norf<err){
00417
00418 lambda+=dlambda;
00419 print_stept(lcid, i+1, lambda, f);
00420 print_flusht();
00421
00422 if (lambda>1.00)
00423 break;
00424
00425 modif++;
00426 if (modif>1){
00427 dlambda*=2.0;
00428 if (dlambda>dlambdamax)
00429 dlambda=dlambdamax;
00430
00431 if (Mesprt==1){
00432 fprintf (stdout,"\n increment must be modified (dlambda increase) dlambda=%e",dlambda);
00433 }
00434 modif=0;
00435 }
00436
00437 continue;
00438 }
00439
00440
00441 for (j=0;j<ini;j++){
00442
00443
00444 Tp->ssle->solve_system (Gtt,Kmat,dd,fb,Outt);
00445
00446 for (k=0;k<n;k++){
00447 d[k]+=dd[k];
00448 }
00449
00450
00451 conductivity_matrix (lcid);
00452 Kmat->gmxv (d,fb);
00453
00454 solution_correction ();
00455 compute_req_valt (0);
00456 approximation ();
00457
00458
00459
00460
00461
00462 for (k=0;k<n;k++){
00463 fb[k]=f[k]-fb[k];
00464 }
00465
00466
00467 norf=normv(fb,n);
00468
00469
00470 fprintf (stdout,"\n increment %ld inner loop j %ld norf=%e dlambda %e",i,j,norf,dlambda);
00471
00472
00473 if (norf<err) break;
00474 }
00475
00476 modif=0;
00477
00478 if (j==ini || norf>err){
00479 dlambda/=2.0;
00480 if (dlambda<dlambdamin){
00481 dlambda=dlambdamin;
00482 break;
00483 }
00484
00485 if (Mesprt==1){
00486 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
00487
00488 }
00489
00490 for (j=0;j<n;j++){
00491 d[j]=db[j];
00492 }
00493 lambda=blambda;
00494 }
00495 else{
00496 lambda+=dlambda;
00497 print_stept(lcid, i+1, lambda, f);
00498 print_flusht();
00499
00500 if (lambda>1.0)
00501 break;
00502 }
00503
00504 }
00505
00506 print_closet();
00507 }