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