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