00001 #include "pnewtonraph.h"
00002 #include "pmtsolver.h"
00003 #include "pglobal.h"
00004 #include "seqfilesm.h"
00005 #include "genfile.h"
00006 #include <string.h>
00007 #include "mpi.h"
00008
00009
00010 void par_solve_timemech_prob ()
00011 {
00012 long lcid;
00013
00014
00015
00016 lcid=0;
00017
00018
00019
00020 par_solve_timemech_prob2 (lcid);
00021
00022 }
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 void par_solve_timemech_prob2 (long lcid)
00035 {
00036 long i, li, nsts, ret;
00037 double dt, dtmin, dtmax,dtdef, end_time;
00038 mt_glob_vec mt_gv;
00039
00040
00041
00042
00043 par_visco_solver_init(lcid, mt_gv);
00044
00045
00046 dt = Mp->timecon.initialtimeincr ();
00047
00048 dtmin=Mp->timecon.dtmin;
00049
00050 dtmax=Mp->timecon.dtmax;
00051
00052 end_time = Mp->timecon.endtime ();
00053
00054
00055 li = i = mt_gv.istep;
00056
00057 nsts=0;
00058
00059 ret = 0;
00060
00061
00062
00063
00064 do
00065 {
00066
00067 i++;
00068
00069
00070 if((Mp->timecon.timefun.tfunc != stat) && Mp->timecon.tct == 0)
00071 dt = Mp->timecon.actualbacktimeincr ();
00072
00073
00074 ret = par_one_step(lcid, Mp->timecon.newtime(dt), dt, i, li, mt_gv);
00075
00076 if (ret >= 0)
00077 {
00078 if (ret == 0)
00079 nsts++;
00080 dtdef = Mp->timecon.actualforwtimeincr();
00081 if (nsts==2)
00082 {
00083 dt*=2.0;
00084 nsts=0;
00085 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00086 }
00087 if (dt>dtdef)
00088 dt=dtdef;
00089
00090 if (dt>dtmax)
00091 dt = dtmax;
00092 }
00093 else
00094 {
00095
00096
00097 dt/=2.0;
00098 Mp->timecon.oldtime ();
00099 i--;
00100
00101 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00102 if (dt<dtmin)
00103 {
00104 if (Myrank==0 && Mespr==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00105 if (Myrank==0 && Mespr==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00106 if (Myrank==0 && Mespr==1) fprintf (stderr,"\n FORCED output of results from this step is performed\n");
00107 compute_req_val (lcid);
00108 print_step_forced(lcid, i, Mp->time, mt_gv.fb);
00109
00110 print_flush();
00111 break;
00112 }
00113 }
00114 }while(Mp->time<=end_time);
00115
00116 print_close ();
00117 }
00118
00119
00120
00121
00122
00123
00124
00125 void par_visco_solver_init(long lcid, mt_glob_vec &mt_gv)
00126 {
00127 long n;
00128 double dt;
00129
00130 n=Ndofm;
00131
00132 mt_gv.alloc(n);
00133
00134 mt_gv.r = Lsrs->give_lhs (lcid);
00135
00136 mt_gv.f = Lsrs->give_rhs (lcid);
00137
00138
00139 Mp->time=Mp->timecon.starttime ();
00140
00141
00142 if (Mp->hdbcont.restore_stat())
00143 {
00144 if (Mespr==1)
00145 fprintf (stdout,"\n Reading of backup file\n");
00146 solver_restore (mt_gv.r, mt_gv.fp, mt_gv.istep, Mp->time, dt, &Mp->timecon, n);
00147
00148 Mm->initmaterialmodels(lcid);
00149 print_init(-1, "at",Pmp->fni,Pmp->fei);
00150
00151 print_flush();
00152 }
00153 else
00154 {
00155 mt_gv.istep=0;
00156
00157 Mm->initmaterialmodels(lcid);
00158
00159 if (Mp->eigstrains > 0)
00160 internal_forces(lcid, mt_gv.fp);
00161
00162 print_init(-1, "wt",Pmp->fni,Pmp->fei);
00163 print_step(lcid, mt_gv.istep, Mp->time, mt_gv.f);
00164
00165 print_flush();
00166 }
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 long par_one_step (long lcid,double time, double dt, long istep, long li, mt_glob_vec &mt_gv)
00192 {
00193 long j,n,ini;
00194 double *f,*fl,*fi,*fb,*fp,*r,*dr,*lhsb;
00195 double norf, err, zero;
00196 matrix lsm_a(3,3);
00197 vector lsm_r(3), lsm_l(3);
00198 Mp->nlman->divc_step = new double[10];
00199 Mp->nlman->divc_err = new double[10];
00200
00201
00202 Mp->time = time;
00203
00204 Mp->istep = istep;
00205
00206
00207 Mp->strainstate=0;
00208
00209
00210 Mp->stressstate=0;
00211
00212
00213 Mp->otherstate=0;
00214
00215
00216 r = mt_gv.r;
00217
00218 f = mt_gv.f;
00219
00220 dr = mt_gv.dr;
00221
00222 fp = mt_gv.fp;
00223
00224 fl = mt_gv.fl;
00225
00226 fi = mt_gv.fi;
00227
00228 fb = mt_gv.fb;
00229
00230 lhsb = mt_gv.lhsb;
00231
00232
00233 n=Ndofm;
00234
00235
00236 ini = Mp->niilnr;
00237
00238
00239 err = Mp->errnr;
00240
00241
00242 zero = Mp->zero;
00243
00244
00245 for (j=0;j<n;j++)
00246 lhsb[j]=r[j];
00247
00248
00249 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n\n -------------------------------------------------------------");
00250 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n MEFEL Time step = %ld, Time %e, Time increment = %e",istep,Mp->time,dt);
00251 if (Myrank==0 && Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------\n");
00252
00253
00254 mefel_right_hand_side (lcid,f,fl);
00255
00256
00257
00258 Mb->comp_sum (f);
00259
00260
00261 fflush(Out);
00262
00263
00264 incr_internal_forces (lcid,fi);
00265
00266
00267
00268 subv(f, fp, fb, n);
00269 addv(fb, fi, fb, n);
00270
00271 norf = par_gnewton_raphson_one_step(lcid, Mp->nlman, fl, r, fb, dr, fi, istep-1, j, li, ini, err);
00272
00273 if (norf>err)
00274 {
00275
00276 if (Mespr==1) fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00277
00278
00279
00280 copyv(lhsb, r, n);
00281
00282 return -1;
00283 }
00284
00285
00286 if ((Mespr==1) && (j > 0)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00287
00288
00289
00290 copyv(f, fp, n);
00291
00292
00293 Mb->comp_sum (fi);
00294
00295
00296 fflush(Out);
00297
00298
00299 Mm->updateipval();
00300 compute_req_val (lcid);
00301 print_step(lcid, istep, Mp->time, fl);
00302
00303 print_flush();
00304
00305 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.save_stat()))
00306 {
00307 if (Mespr==1)
00308 fprintf (stdout,"\n Creating backup file\n");
00309 solver_save (r,fp,istep,Mp->time,dt,&Mp->timecon,n);
00310 }
00311 if (j == 0)
00312 return 0;
00313
00314 return 1;
00315 }
00316
00317
00318
00319
00320
00321
00322
00323
00324 void par_solve_timemech_prob (long lcid)
00325 {
00326 long i,j,k,n,ini,nsts,nii;
00327 double end_time,zero,err,dt,dtmin,dtdef,norfa,norf;
00328 double *f,*fi,*fp,*fb,*fl,*dr,*r,*lhsb;
00329 matrix lsm_a(3,3);
00330 vector lsm_r(3), lsm_l(3);
00331 double lsm_res;
00332
00333
00334
00335
00336 ini = Mp->niilnr;
00337
00338 err = Mp->errnr;
00339
00340 n = Ndofm;
00341
00342 zero = Mp->zero;
00343
00344
00345
00346 r = Lsrs->give_lhs (lcid);
00347
00348 f = Lsrs->give_rhs (lcid);
00349
00350
00351
00352
00353 dr = new double [n];
00354 nullv (dr,n);
00355
00356 fp = new double [n];
00357 nullv (fp,n);
00358
00359 fl = new double [n];
00360 nullv (fl,n);
00361
00362 fi = new double [n];
00363 nullv (fi,n);
00364
00365 fb = new double [n];
00366 nullv (fb,n);
00367
00368 lhsb = new double [n];
00369 nullv (lhsb,n);
00370
00371
00372 Mp->time=Mp->timecon.starttime ();
00373
00374 end_time = Mp->timecon.endtime ();
00375
00376 dt = Mp->timecon.initialtimeincr ();
00377
00378 dtmin=Mp->timecon.dtmin;
00379
00380
00381 i=0;
00382
00383 nsts=0;
00384
00385
00386
00387 if (Mp->hdbcont.restore_stat()){
00388 if (Mespr==1)
00389 fprintf (stdout,"\n Reading of backup file\n");
00390 solver_restore (r,fp,i,Mp->time,dt,&Mp->timecon,n);
00391
00392 Mm->initmaterialmodels(lcid);
00393 print_init(-1, "at",Pmp->fni,Pmp->fei);
00394 print_step(lcid, i, Mp->time, f);
00395 print_flush();
00396 }
00397 else
00398 {
00399
00400 Mm->initmaterialmodels(lcid);
00401 print_init(-1, "wt",Pmp->fni,Pmp->fei);
00402 print_step(lcid, i, Mp->time, f);
00403 print_flush();
00404 }
00405
00406
00407
00408
00409
00410
00411 do{
00412
00413 Mp->time = Mp->timecon.newtime (dt);
00414
00415 i++;
00416
00417
00418 Mp->strainstate=0;
00419
00420
00421 Mp->stressstate=0;
00422
00423
00424 Mp->otherstate=0;
00425
00426
00427 for (j=0;j<n;j++){
00428 lhsb[j]=r[j];
00429 }
00430
00431 if ((Myrank==0)&&(Mespr==1)) fprintf (stdout,"\n\n -------------------------------------------------------------");
00432 if ((Myrank==0)&&(Mespr==1)) fprintf (stdout,"\n Time step = %ld, Time %e, Time increment = %e",i,Mp->time,dt);
00433 if ((Myrank==0)&&(Mespr==1)) fprintf (stdout,"\n -------------------------------------------------------------\n");
00434
00435
00436 mefel_right_hand_side (lcid,f,fl);
00437
00438
00439
00440 Mb->comp_sum (f);
00441
00442
00443 fflush(Out);
00444
00445
00446 incr_internal_forces (lcid,fi);
00447
00448
00449 for (j=0;j<n;j++){
00450 fb[j]=f[j]-fp[j]+fi[j];
00451 }
00452
00453
00454 if ((Mp->nlman->stmat==tangent_stiff) || (i == 1) || (Smat==NULL))
00455 stiffness_matrix (lcid);
00456
00457
00458 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00459
00460
00461 for (j=0;j<n;j++){
00462 r[j]+=dr[j];
00463 }
00464
00465
00466 internal_forces (lcid,fi);
00467
00468
00469 norfa=Psolm->pss (f,f,Out);
00470
00471 for (j=0;j<n;j++){
00472 fb[j] = fl[j] - fi[j];
00473 }
00474
00475 norf=Psolm->pss(fb,fb,Out);
00476
00477 if (norfa > 0.0)
00478 norf /= norfa;
00479
00480 if ((Myrank==0)&&(Mespr==1))
00481 fprintf (stdout,"\n\n Norf before inner iteration loop norf = %e\n\n\n",norf);
00482 j = 0;
00483 if (Psolm->compare_on_master(norf, err) <= 0)
00484 {
00485
00486 nsts++;
00487
00488 nii=1;
00489 }
00490 else
00491 nii=0;
00492
00493
00494 if (nii==0){
00495
00496
00497 fillm(0.0, lsm_a);
00498 fillv(0.0, lsm_r);
00499 for (j=0;j<ini;j++)
00500 {
00501
00502 if ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0))
00503 stiffness_matrix (lcid);
00504
00505 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00506
00507
00508 for (k=0;k<n;k++)
00509 r[k]+=dr[k];
00510
00511
00512 internal_forces (lcid,fi);
00513
00514
00515 for (k=0;k<n;k++)
00516 fb[k]=fl[k]-fi[k];
00517
00518
00519 norf=Psolm->pss(fb,fb,Out);
00520 if (norfa > 0.0)
00521 norf /= norfa;
00522
00523
00524 if ((Myrank==0) && (Mespr==1))
00525 fprintf (stdout,"\n Time increment %lf inner loop j %ld norf=%e",Mp->time,j,norf);
00526 if (Psolm->compare_on_master(norf, err) <= 0)
00527 {
00528
00529
00530 nsts++;
00531
00532 nii=1;
00533 break;
00534 }
00535 else
00536 nii=0;
00537
00538
00539 if (j > 1)
00540 {
00541 lsm_res = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,1);
00542 if (Psolm->compare_on_master(lsm_res, 0.0) > 0)
00543 {
00544 if ((Myrank==0) && (Mespr==1))
00545 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00546
00547 nii=0;
00548 break;
00549 }
00550 }
00551 else
00552 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
00553 }
00554 }
00555
00556 if (nii==0)
00557 {
00558
00559 if ((Myrank==0) && (Mespr==1))
00560 fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00561
00562
00563 for (j=0;j<n;j++){
00564 r[j]=lhsb[j];
00565 }
00566
00567
00568
00569 dt/=2.0;
00570 Mp->timecon.oldtime ();
00571 i--;
00572
00573 if ((Myrank==0) && (Mespr==1))
00574 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00575 if (Psolm->compare_on_master(dt, dtmin) < 0){
00576 fflush(stdout);
00577 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00578 if ((Myrank==0) && (Mespr==1)) fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
00579 compute_req_val (lcid);
00580 print_step_forced(lcid, i, Mp->time, f);
00581 print_flush();
00582 break;
00583 }
00584 }
00585 else{
00586
00587 if ((Myrank==0) && (Mespr==1)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00588
00589
00590 for (j=0;j<n;j++){
00591 fp[j]=f[j];
00592 }
00593
00594 dtdef = Mp->timecon.actualforwtimeincr ();
00595 if (nsts==2){
00596 dt*=2.0;
00597 nsts=0;
00598 if ((Myrank==0) && (Mespr==1))
00599 fprintf (stdout,"\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00600 }
00601 if (dt>dtdef)
00602 dt=dtdef;
00603
00604
00605 Mb->comp_sum (fi);
00606 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fix=%le, fiy=%le, fiz=%le",
00607 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00608 fflush(Out);
00609
00610
00611 Mm->updateipval();
00612 compute_req_val (lcid);
00613 print_step(lcid, i, Mp->time, f);
00614 print_flush();
00615
00616 if ((Mp->timecon.isitimptime ()==1) && (Mp->hdbcont.hdbtype))
00617 {
00618 if ((Myrank==0) && (Mespr==1))
00619 fprintf (stdout,"\n Creating backup file\n");
00620 solver_save (r,fp,i,Mp->time,dt,&Mp->timecon,n);
00621 }
00622 }
00623
00624 }while(Mp->time<=end_time);
00625
00626 print_close ();
00627
00628
00629 delete [] fi;
00630 delete [] fb;
00631 delete [] fp;
00632 delete [] dr;
00633 delete [] fl;
00634 delete [] lhsb;
00635
00636
00637 }