00001 #include "pmtsolver.h"
00002 #include "pglobal.h"
00003 #include "seqfilesm.h"
00004 #include "genfile.h"
00005 #include <string.h>
00006 #include "mpi.h"
00007
00008
00009
00010
00011 void par_solve_timemech_prob ()
00012 {
00013 long lcid;
00014
00015
00016
00017 lcid=0;
00018
00019
00020 par_solve_timemech_prob (lcid);
00021
00022 }
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 void par_solve_timemech_prob (long lcid)
00033 {
00034 long i,j,k,n,ini,nsts,nii;
00035 double end_time,zero,err,dt,dtmin,dtdef,norfa,norf;
00036 double *f,*fi,*fp,*fb,*fl,*dr,*r,*lhsb;
00037 matrix lsm_a(3,3);
00038 vector lsm_r(3), lsm_l(3);
00039 double lsm_res;
00040
00041
00042
00043 ini = Mp->niilnr;
00044
00045 err = Mp->errnr;
00046
00047 n = Ndofm;
00048
00049 zero = Mp->zero;
00050
00051
00052 fp = new double [n];
00053 fb = new double [n];
00054 fl = new double [n];
00055 fi = new double [n];
00056 dr = new double [n];
00057 lhsb = new double [n];
00058
00059
00060 nullv (fp,n);
00061 nullv (fb,n);
00062 nullv (fl,n);
00063 nullv (fi,n);
00064 nullv (dr,n);
00065 nullv (lhsb,n);
00066
00067
00068 r = Lsrs->give_lhs (lcid);
00069
00070 f= Lsrs->give_rhs (lcid);
00071
00072
00073 Mp->time=Mp->timecon.starttime ();
00074
00075 end_time = Mp->timecon.endtime ();
00076
00077 dt = Mp->timecon.initialtimeincr ();
00078
00079 dtmin=Mp->timecon.dtmin;
00080
00081
00082
00083 Gtm->comp_domain_sizes ();
00084 Mb->alloc_sumcomp ();
00085
00086
00087
00088 nullv (fp,n);
00089
00090 i=0;
00091
00092 nsts=0;
00093
00094
00095
00096 if (Mp->hdbcont.restore_stat()){
00097 solver_restore (r,fp,i,Mp->time,dt,&Mp->timecon,n);
00098
00099 Mm->initmaterialmodels();
00100 print_init(-1, "at",Pmp->fni,Pmp->fei);
00101 print_step(lcid, i, Mp->time, f);
00102 }
00103 else
00104 {
00105
00106 Mm->initmaterialmodels();
00107 print_init(-1, "wt",Pmp->fni,Pmp->fei);
00108 print_step(lcid, i, Mp->time, f);
00109 print_flush();
00110 }
00111
00112
00113
00114
00115
00116 do{
00117
00118
00119 Mp->time = Mp->timecon.newtime (dt);
00120
00121 i++;
00122
00123
00124 Mp->strainstate=0;
00125
00126
00127 Mp->stressstate=0;
00128
00129
00130 Mp->otherstate=0;
00131
00132
00133
00134 for (j=0;j<n;j++){
00135 lhsb[j]=r[j];
00136 }
00137
00138 if (Mespr==1)
00139 fprintf (stdout,"\n time %e",Mp->time);
00140 fprintf (Out,"\n time %e",Mp->time);
00141
00142
00143 mefel_right_hand_side (lcid,f,fl);
00144
00145 Mb->comp_sum (f);
00146 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fx=%le, fy=%le, fz=%le",
00147 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00148 fflush(Out);
00149
00150
00151 incr_internal_forces (lcid,fi);
00152
00153
00154 for (j=0;j<n;j++){
00155 fb[j]=f[j]-fp[j]+fi[j];
00156 }
00157
00158
00159 if ((Mp->stmat==tangent_stiff) || (i == 1) || (Smat==NULL))
00160 stiffness_matrix (lcid);
00161
00162
00163
00164 Psol->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00165
00166
00167 for (j=0;j<n;j++){
00168 r[j]+=dr[j];
00169 }
00170
00171
00172 internal_forces (lcid,fi);
00173
00174 norfa=Psol->pss (f,Out);
00175 for (j=0;j<n;j++){
00176 fb[j] = fl[j] - fi[j];
00177 }
00178
00179 norf=Psol->pss(fb,Out);
00180
00181 j = 0;
00182 if (norfa > 0.0)
00183 norf /= norfa;
00184 if (Psol->compare_on_master(norf, err) <= 0)
00185 {
00186
00187 nsts++;
00188
00189 nii=1;
00190 }
00191 else
00192 nii=0;
00193
00194
00195
00196 if (nii==0){
00197
00198
00199 fillm(0.0, lsm_a);
00200 fillv(0.0, lsm_r);
00201 for (j=0;j<ini;j++){
00202
00203
00204
00205
00206 Psol->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr);
00207 for (k=0;k<n;k++)
00208 r[k]+=dr[k];
00209
00210 internal_forces (lcid,fi);
00211
00212 for (k=0;k<n;k++)
00213 fb[k]=fl[k]-fi[k];
00214
00215 norf=Psol->pss(fb,Out);
00216 if (norfa > 0.0)
00217 norf /= norfa;
00218 if ((Myrank==0) && (Mespr==1))
00219 fprintf (stdout,"\n Time increment %lf inner loop j %ld norf=%e",Mp->time,j,norf);
00220 if (Psol->compare_on_master(norf, err) <= 0)
00221 {
00222
00223 nsts++;
00224
00225 nii=1;
00226 break;
00227 }
00228 else
00229 nii=0;
00230
00231
00232 if (j > 1)
00233 {
00234 lsm_res = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, Mp->zero,1);
00235 if (Psol->compare_on_master(lsm_res, 0.0) > 0)
00236 {
00237 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00238
00239 nii=0;
00240 break;
00241 }
00242 }
00243 else
00244 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, Mp->zero,0);
00245 }
00246 }
00247
00248 if (nii==0){
00249
00250 for (j=0;j<n;j++){
00251 r[j]=lhsb[j];
00252 }
00253
00254
00255 dt/=2.0;
00256 Mp->timecon.oldtime ();
00257 i--;
00258
00259 if ((Myrank==0) && (Mespr==1))
00260 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00261 if (dt<dtmin){
00262 fprintf (stderr,"\n\n time increment is less than minimum time increment");
00263 fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
00264 compute_req_val(lcid);
00265 print_step_forced(lcid, i, Mp->time, f);
00266 print_flush();
00267 abort();
00268 }
00269 }
00270 else{
00271
00272 if ((Myrank==0) && (Mespr==1))
00273 fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00274
00275 for (j=0;j<nm;j++){
00276 fp[j]=f[j];
00277 }
00278 dtdef = Mp->timecon.actualforwtimeincr ();
00279 if (nsts==2){
00280 dt*=2.0;
00281 nsts=0;
00282 if ((Myrank==0) && (Mespr==1))
00283 fprintf (stdout,"\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00284 }
00285 if (dt>dtdef)
00286 dt=dtdef;
00287
00288
00289 Mb->comp_sum (fi);
00290 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fix=%le, fiy=%le, fiz=%le",
00291 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00292 Mm->updateipval();
00293
00294 compute_req_val (lcid);
00295 print_step(lcid, i, Mp->time, f);
00296 print_flush();
00297 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.hdbtype)){
00298 if ((Myrank==0) && (Mespr==1))
00299 fprintf (stdout,"\n Creating backup file for MEFEL\n");
00300 solver_save (r,f,ii,Mp->time,dt,&Mp->timecon,Ndofm);
00301 }
00302 }
00303 }while(Mp->time<end_time);
00304
00305 print_close ();
00306
00307 delete [] dr;
00308 delete [] fi;
00309 delete [] fb;
00310 delete [] fp;
00311 delete [] fl;
00312 delete [] lhsb;
00313 }
00314
00315
00316
00317 void par_solve_timemech_prob_old ()
00318 {
00319 long i,j,n,lcid;
00320 double dt,end_time,zero,*f,*fi,*fp,*dr,*lhs,*rhs;
00321 char file[FNAMELEN+20];
00322 FILE *aux;
00323 char *path;
00324 char *name;
00325 char *suffix;
00326
00327
00328 Mp->filename_decomposition (Outdm->outfn,path,name,suffix);
00329 sprintf (file,"%s.bac",Outdm->outfn);
00330 aux = fopen (file,"w");
00331
00332
00333 lcid=0;
00334
00335
00336 n = Ndofm;
00337
00338 zero = Mp->zero;
00339
00340 fp = new double [n];
00341 f = new double [n];
00342 fi = new double [n];
00343 dr = new double [n];
00344 memset (fp,0,n*sizeof(double));
00345 memset (f,0,n*sizeof(double));
00346 memset (fi,0,n*sizeof(double));
00347 memset (dr,0,n*sizeof(double));
00348
00349
00350 lhs = Lsrs->give_lhs (lcid);
00351 rhs = Lsrs->give_rhs (lcid);
00352
00353
00354
00355 Mp->time=Mp->timecon.starttime ();
00356
00357 end_time = Mp->timecon.endtime ();
00358
00359
00360
00361
00362
00363 Gtm->comp_domain_sizes ();
00364 Mb->alloc_sumcomp ();
00365
00366
00367
00368 nullv (fp,n);
00369 i=0;
00370
00371
00372 print_init(-1, "wt",Pmp->fni,Pmp->fei);
00373 print_step(lcid, i, Mp->time, f);
00374 print_flush();
00375
00376
00377
00378
00379
00380
00381
00382 do{
00383
00384
00385 Mm->updateipval();
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 mefel_right_hand_side (lcid,f);
00402
00403
00404
00405
00406
00407
00408 Mb->comp_sum (f);
00409
00410
00411
00412 Mp->phase=1;
00413 internal_forces (lcid,fi);
00414
00415
00416 for (j=0;j<n;j++){
00417 rhs[j]=f[j]-fp[j]+fi[j];
00418 fp[j]=f[j];
00419 }
00420
00421
00422
00423 stiffness_matrix (lcid);
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434 Psol->par_linear_solver (Gtm,Smat,dr,rhs,Out,Mespr);
00435
00436
00437
00438
00439
00440
00441 for (j=0;j<n;j++){
00442 lhs[j]+=dr[j];
00443 }
00444
00445
00446 Mp->phase=2;
00447 internal_forces (lcid,fi);
00448
00449 compute_req_val (lcid);
00450
00451
00452
00453 Mp->time = Mp->timecon.newtime ();
00454 i++;
00455
00456
00457 if (Myrank==0)
00458 fprintf (stdout,"\n step %ld time %lf",i,Mp->time);
00459
00460
00461 if (Mp->timecon.isitimptime ()==1){
00462 solver_save (lhs,fp,i,Mp->time,dt,&Mp->timecon,n);
00463 }
00464
00465
00466
00467
00468 print_step(lcid, i, Mp->time, f);
00469 print_flush();
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 }while(Mp->time<end_time);
00491
00492
00493 print_close ();
00494
00495 fclose (aux);
00496
00497 delete [] dr; delete [] fi; delete [] fp; delete [] f;
00498
00499 }
00500