00001 #include "ppcsolver.h" 00002 00003 #include "global.h" 00004 #include "globalt.h" 00005 #include "globalc.h" 00006 #include "pglobalc.h" 00007 00008 #include "globmat.h" 00009 #include "globmatt.h" 00010 #include "globmatc.h" 00011 00012 #include "mechprint.h" 00013 #include "transprint.h" 00014 00015 /** 00016 function solves partially coupled thermo-hydro-mechanical problem in parallel 00017 transport problems are linear or nonlinear 00018 mechanical analysis is always nonlinear 00019 00020 JK 00021 */ 00022 void par_solve_pcouplprob () 00023 { 00024 long lcid; 00025 00026 // load case id must be equal to zero 00027 // the problems are nonlinear and superposition method cannot be used 00028 lcid=0; 00029 00030 switch (Cp->tnlinsol){ 00031 case newtonc:{ 00032 par_newton_raphson_parcoupl (lcid); 00033 break; 00034 } 00035 default:{ 00036 print_err("unknown solver of nonlinear equation system is required",__FILE__,__LINE__,__func__); 00037 } 00038 } 00039 00040 } 00041 00042 00043 /** 00044 function solves partially coupled thermo-hydro-mechanical time-dependent problem 00045 by the Newton-Raphson method 00046 00047 @param lcid - load case id 00048 00049 JK, 5.8.2008 00050 */ 00051 void par_newton_raphson_parcoupl (long lcid) 00052 { 00053 switch (Tp->tprob){ 00054 case nonstationary_problem:{ 00055 // linear nonstationary transport problem 00056 // nonlinear mechanical problem 00057 par_newton_raphson_parcoupl_lin_vform (lcid); 00058 break; 00059 } 00060 case nonlinear_nonstationary_problem:{ 00061 // nonlinear nonstationary transport problem 00062 // nonlinear mechanical problem 00063 par_newton_raphson_parcoupl_nonlin (lcid); 00064 break; 00065 } 00066 default:{ 00067 print_err("unknown TRFEL problem type is required",__FILE__,__LINE__,__func__); 00068 } 00069 } 00070 } 00071 00072 /** 00073 function solves partially coupled problems in parallel, 00074 transport processes are linear, mechanical analysis is nonlinear 00075 systems of nonlinear equations are solved by the Newton-Raphson method 00076 trapezoidal rule is in the v-form, it means, nodal time derivatives are calculated 00077 nodal values are computed later from time derivatives 00078 00079 @param lcid - load case id 00080 00081 JK, 19.5.2004, 11.6.2008 00082 */ 00083 void par_newton_raphson_parcoupl_lin_vform (long lcid) 00084 { 00085 long i,j,k,nm,nt,ini,nii,nsts; 00086 double zero,dt,dtmin,dtdef,end_time,alpha,norfa,norf,err,s; 00087 double *r,*dr,*d,*f,*fb,*fp,*fi,*fl,*p,*lhsb,*lhst,*lhstb,*tdlhst,*rhst; 00088 matrix lsm_a(3,3); 00089 vector lsm_r(3), lsm_l(3); 00090 MPI_Status stat; 00091 char file[FNAMELEN+20]; 00092 FILE *aux; 00093 char *path; 00094 char *name; 00095 char *suffix; 00096 00097 00098 // number of mechanical degrees of freedom 00099 nm=Ndofm; 00100 // number of transport degrees of freedom 00101 nt=Ndoft; 00102 00103 //Mp->filename_decomposition (Outdm->outfn,path,name,suffix); 00104 //sprintf (file,"%s.bac",Outdm->outfn); 00105 //aux = fopen ("zaloha","w"); 00106 00107 00108 // *************************** 00109 // vectors of transport part 00110 // *************************** 00111 00112 // vector of nodal values 00113 lhst=Lsrst->give_lhs (lcid); 00114 // vector of time derivatives of nodal values 00115 tdlhst=Lsrst->give_tdlhs (lcid); 00116 // vector of the right hand side 00117 rhst=Lsrst->give_rhs (lcid); 00118 00119 // array containing predictors 00120 d = new double [nt]; 00121 nullv (d,nt); 00122 // predictor vector 00123 p = new double [nt]; 00124 nullv (p,nt); 00125 // backup of the lhs vector 00126 lhstb = new double [nt]; 00127 nullv (lhstb,nt); 00128 00129 00130 // **************************** 00131 // vectors of mechanical part 00132 // **************************** 00133 00134 // vector of nodal displacements 00135 r = Lsrs->give_lhs (lcid); 00136 // vector of nodal prescribed forces 00137 f = Lsrs->give_rhs (lcid); 00138 00139 // vector of increments of nodal displacements 00140 dr = new double [nm]; 00141 nullv (dr,nm); 00142 // vector of prescribed forces from the previous step 00143 fp = new double [nm]; 00144 nullv (fp,nm); 00145 // vector of force loads 00146 fl = new double [nm]; 00147 nullv (fl,nm); 00148 // vector of internal forces 00149 fi = new double [nm]; 00150 nullv (fi,nm); 00151 // auxiliary force vector 00152 fb = new double [nm]; 00153 nullv (fb,nm); 00154 // backup of the nodal displacements 00155 lhsb = new double [nm]; 00156 nullv (lhsb,nm); 00157 00158 00159 // coefficient of trapezoidal method 00160 alpha=Tp->alpha; 00161 // computer zero 00162 zero=Cp->zero; 00163 // maximum number of iterations in inner loop in mechanics 00164 ini = Mp->niilnr; 00165 // required norm of vector of unbalanced forces 00166 err = Mp->errnr; 00167 00168 00169 00170 // starting time 00171 Cp->time=Cp->timecon.starttime (); 00172 // time increment 00173 dt=Cp->timecon.initialtimeincr (); 00174 // end time 00175 end_time = Cp->timecon.endtime (); 00176 00177 00178 // mechanical time controller is rewritten by coupled time controller 00179 Mp->timecon.take_values (Cp->timecon); 00180 //Mp->timecon.seconds_days (); 00181 Mp->time=Cp->time; 00182 00183 // transport time controller is rewritten by coupled time controller 00184 Tp->timecont.take_values (Cp->timecon); 00185 Tp->time = Cp->time; 00186 00187 00188 // nodes - integration points interpolation 00189 approximation (); 00190 // initial temperature is approximated to the mechanical integration points 00191 approximation_inittemper (); 00192 approximation_humid (); 00193 00194 actual_previous_nodval (); 00195 00196 Gtm->comp_domain_sizes (); 00197 Mb->alloc_sumcomp (); 00198 00199 00200 // initiation of mechanical material models 00201 Mm->initmaterialmodels(); 00202 // initiation of transport material models 00203 Tm->initmaterialmodels(); 00204 00205 compute_req_valt (0); 00206 00207 00208 // number of time step 00209 i=0; 00210 // number of successful time steps (number of steps without inner iterations) 00211 nsts=0; 00212 00213 // **************************************** 00214 // print of initial state of the problem 00215 // **************************************** 00216 print_init(-1, "wt",Pcp->fnim,Pcp->feim); 00217 print_step(0,i,Cp->time,NULL); 00218 print_flush(); 00219 00220 print_initt(-1, "wt",Pcp->fnit,Pcp->feit); 00221 print_stept(0,i,Cp->time,NULL); 00222 print_flusht(); 00223 00224 00225 // ************************** 00226 // main iteration loop *** 00227 // ************************** 00228 00229 do{ 00230 00231 // update of time 00232 Cp->time=Cp->timecon.newtime (dt); 00233 00234 if (Mesprc != 0) 00235 fprintf (stdout,"\n\n time %le, step number %ld",Cp->time,i); 00236 00237 // ******************************** 00238 // transport part of the problem 00239 // ******************************** 00240 00241 00242 // conductivity matrix 00243 conductivity_matrix (lcid); 00244 00245 // capacity matrix 00246 capacity_matrix (lcid); 00247 00248 // predictor d+(1-alpha)*dt*v 00249 for (j=0;j<nt;j++){ 00250 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j]; 00251 } 00252 00253 // right hand side of transport part 00254 trfel_right_hand_side (0,rhst,nt); 00255 00256 // auxiliary vector 00257 Kmat->gmxv (d,p); 00258 00259 // matrix of the system of equations 00260 // C + alpha.dt.K 00261 Kmat->scalgm (dt*alpha); 00262 Kmat->addgm (1.0,*Cmat); 00263 00264 for (j=0;j<nt;j++){ 00265 rhst[j] = rhst[j] - d[j]; 00266 d[j]=tdlhst[j]; 00267 } 00268 00269 // solution of the system of algebraic equations 00270 Psolt->par_linear_solver (Gtt,Kmat,tdlhst,rhst,Outt,Mespr); 00271 00272 // backup of attained nodal values 00273 for (j=0;j<nt;j++){ 00274 lhstb[j]=lhst[j]; 00275 } 00276 00277 // new nodal values 00278 for (j=0;j<nt;j++){ 00279 s=(1.0-alpha)*d[j]+alpha*tdlhst[j]; 00280 lhst[j]+=dt*s; 00281 } 00282 00283 // physically corrected solution 00284 solution_correction (); 00285 // approximation of nodal values into ontegration points 00286 approximation (); 00287 00288 actual_previous_nodval (); 00289 00290 00291 00292 print_stept(0,i,Cp->time,NULL); 00293 print_flusht(); 00294 00295 00296 // ********************************************************* 00297 // end of transport part and beginning of mechanical part 00298 // ********************************************************* 00299 00300 if (Mp->time <= Tp->time){ 00301 00302 // indicator of strain computation 00303 // it means, strains at integration points have not been computed 00304 Mp->strainstate=0; 00305 // indicator of stress computation 00306 // it means, stresses at integration points have not been computed 00307 Mp->stressstate=0; 00308 // indicator of computation of other array 00309 // it means, stresses at integration points have not been computed 00310 Mp->otherstate=0; 00311 00312 // approximation of temperatures from TRFEL nodes to MEFEL integration points 00313 approximation_temper (); 00314 // approximation of humidities from TRFEL nodes to MEFEL integration points 00315 approximation_humid (); 00316 00317 00318 // backup of attained nodal values 00319 for (j=0;j<nm;j++){ 00320 lhsb[j]=r[j]; 00321 } 00322 00323 // assembling of the right hand side (prescribed load) 00324 mefel_right_hand_side (lcid,f,fl); 00325 00326 // computation of sums of force components in particular directions 00327 // it is used in special creep and consolidation models 00328 Mb->comp_sum (f); 00329 00330 // computation of unbalanced forces which will be on the right hand side 00331 incr_internal_forces (lcid,fi); 00332 00333 for (j=0;j<nm;j++){ 00334 fb[j]=f[j]-fp[j]+fi[j]; 00335 fp[j]=f[j]; 00336 } 00337 00338 // assembling of tangent stiffness matrix 00339 stiffness_matrix (lcid); 00340 00341 // solution of K(r).v=F 00342 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr); 00343 00344 // new displacement vector 00345 for (j=0;j<nm;j++){ 00346 r[j]+=dr[j]; 00347 } 00348 00349 // computation of internal forces 00350 internal_forces (lcid,fi); 00351 00352 norfa=Psolm->pss (f,Out); 00353 00354 for (j=0;j<nm;j++){ 00355 fb[j] = fl[j] - fi[j]; 00356 } 00357 // norm of vector of unbalanced forces 00358 norf=Psolm->pss (fb,Out); 00359 00360 if (norfa > 0.0) 00361 norf /= norfa; 00362 00363 if (Psolm->compare_on_master(norf, err) <= 0){ // norf <= err 00364 // number of successful time steps 00365 nsts++; 00366 // no inner iteration 00367 nii=1; 00368 } 00369 else{ 00370 // inner iteration is needed 00371 nii=0; 00372 } 00373 00374 if (nii==0){ 00375 // iteration of unbalanced forces caused by time independent models 00376 // internal iteration loop 00377 for (j=0;j<ini;j++){ 00378 00379 if(Cp->fcsolv == fullnewtonc){ 00380 //cleaning of stiffness matrix 00381 if (Smat != NULL){ 00382 delete Smat; 00383 Smat=NULL; 00384 } 00385 //assembling of tangent stiffness matrix for fullnewton 00386 stiffness_matrix (lcid); 00387 } 00388 00389 // solution of K(r).v=F 00390 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr); 00391 00392 for (k=0;k<nm;k++) 00393 r[k]+=dr[k]; 00394 00395 // computation of internal forces 00396 internal_forces (lcid,fi); 00397 00398 // vector of unbalanced forces 00399 for (k=0;k<nm;k++) 00400 fb[k]=fl[k]-fi[k]; 00401 00402 // norm of vector of unbalanced forces 00403 norf=Psolm->pss(fb,Out); 00404 if (norfa > 0.0) 00405 norf /= norfa; 00406 00407 if ((Myrank==0) && (Mespr==1)) 00408 fprintf (stdout,"\n Time increment %lf inner loop j %ld norf=%e",Mp->time,j,norf); 00409 if (Psolm->compare_on_master(norf, err) <= 0){ // norf <= err 00410 // number of successful time steps 00411 nsts++; 00412 // no inner iteration 00413 nii=1; 00414 } 00415 else 00416 nii=0; 00417 00418 if (nii==1) 00419 break; 00420 } 00421 } 00422 00423 if (nii==0){ 00424 // inner iteration was not successful 00425 // algorithm has to return to the beginning 00426 00427 // backup is used for actual values 00428 for (j=0;j<nm;j++){ 00429 r[j]=lhsb[j]; 00430 } 00431 // reduction of the time increment because 00432 // inner loop was not able to enforce equilibrium 00433 dt/=2.0; 00434 Mp->timecon.oldtime (); 00435 i--; 00436 00437 if ((Myrank==0) && (Mespr==1)) 00438 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium"); 00439 if (dt<dtmin){ 00440 fprintf (stderr,"\n\n time increment is less than minimum time increment"); 00441 fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__); 00442 abort(); 00443 } 00444 00445 } 00446 else{ 00447 // inner iteration was successful 00448 00449 dtdef = Mp->timecon.actualforwtimeincr (); 00450 if (nsts==2){ 00451 dt*=2.0; 00452 nsts=0; 00453 if ((Myrank==0) && (Mespr==1)) 00454 fprintf (stdout,"\n time increment is enlarged because no inner loop was neccessary in previous 3 steps"); 00455 } 00456 if (dt>dtdef) 00457 dt=dtdef; 00458 00459 Mm->updateipval(); 00460 00461 compute_req_val (lcid); 00462 print_step(lcid, i, Mp->time, f); 00463 print_flush(); 00464 } 00465 00466 } 00467 00468 }while(Cp->time<=end_time); 00469 00470 print_close (); 00471 print_closet (); 00472 00473 delete [] fi; 00474 delete [] fb; 00475 delete [] fp; 00476 delete [] f; 00477 delete [] p; 00478 delete [] d; 00479 delete [] dr; 00480 //delete [] r; 00481 00482 //fclose (aux); 00483 00484 } 00485 00486 00487 00488 00489 /** 00490 function solves partially coupled problems in parallel, 00491 transport processes are nonlinear, mechanical analysis is nonlinear 00492 systems of nonlinear equations are solved by the Newton-Raphson method 00493 00494 @param lcid - load case id 00495 00496 JK, 19.9.2007, 11.6.2008 00497 */ 00498 void par_newton_raphson_parcoupl_nonlin (long lcid) 00499 { 00500 /* 00501 long i,ii,j,k,nm,nt,ini; 00502 double zero,dt,end_time,alpha; 00503 double *dr,*fb,*r,*d,*f,*fp,*fi,*p,*lhst,*tdlhst,*rhst; 00504 double *fbt,*fit,*lhst_last; 00505 //double s; 00506 double norf_last; 00507 double norf, norfa, err; 00508 char fname[1020]; 00509 char *path; 00510 char *name; 00511 char *suffix; 00512 matrix lsm_a(3,3); 00513 vector lsm_r(3), lsm_l(3); 00514 FILE *aux; 00515 MPI_Status stat; 00516 00517 00518 Mp->filename_decomposition (Outdm->outfn,path,name,suffix); 00519 Mp->path = path; 00520 Mp->filename=name; 00521 Mp->suffix=suffix; 00522 sprintf (fname,"%s%s.bac",path, name); 00523 aux = fopen (fname,"w"); 00524 00525 00526 // maximum number of iterations in inner loop 00527 ini = Cp->niilnr; 00528 // required norm of vector of unbalanced forces 00529 err = Cp->errnr; 00530 // number of mechanical degrees of freedom 00531 nm=Ndofm; 00532 // number of transport degrees of freedom 00533 nt=Ndoft; 00534 00535 // *************************** 00536 // vectors of transport part 00537 // *************************** 00538 // left hand side 00539 lhst=Lsrst->give_lhs (0); 00540 // time derivatives of nodal values 00541 tdlhst=Lsrst->give_tdlhs (0); 00542 // right hand side 00543 rhst=Lsrst->give_rhs (0); 00544 00545 // array containing predictors 00546 d = new double [nt]; 00547 // auxiliary vector 00548 p = new double [nt]; 00549 fbt = new double [nt]; 00550 fit = new double [nt]; 00551 lhst_last = new double [nt]; 00552 00553 // initial values 00554 nullv (lhst,nt); 00555 nullv (tdlhst,nt); 00556 nullv (d,nt); 00557 nullv (p,nt); 00558 nullv (fbt,nt); 00559 nullv (fit,nt); 00560 00561 // **************************** 00562 // vectors of mechanical part 00563 // **************************** 00564 fp = new double [nm]; 00565 fb = new double [nm]; 00566 fi = new double [nm]; 00567 dr = new double [nm]; 00568 00569 // initial values 00570 nullv (fp,nm); 00571 nullv (fb,nm); 00572 nullv (fi,nm); 00573 nullv (dr,nm); 00574 00575 // initialization phase 00576 r = Lsrs->give_lhs (lcid); 00577 f = Lsrs->give_rhs (lcid); 00578 00579 // coefficient of trapezoidal method 00580 alpha=Tp->alpha; 00581 zero=Tp->zero; 00582 00583 // driving time from METR 00584 // starting time 00585 Cp->time=Cp->timecon.starttime (); 00586 // time increment 00587 dt=Cp->timecon.initialtimeincr (); 00588 // end time 00589 end_time = Cp->timecon.endtime (); 00590 00591 00592 // mechanical time controller is rewritten by coupled time controller 00593 Mp->time = Cp->time; 00594 Mp->timecon.take_values (Cp->timecon); 00595 // transport time controller is rewritten by coupled time controller 00596 Tp->time = Cp->time; 00597 Tp->timecont.take_values (Cp->timecon); 00598 00599 // nodes - integration points interpolation in TRFEL 00600 // values at integration points are interpolated from nodal values 00601 approximation (); 00602 // initial temperature is approximated to the mechanical integration points 00603 approximation_inittemper (); 00604 00605 approximation_humid (); 00606 00607 // initiation of mechanical material models 00608 Mm->initmaterialmodels(); 00609 // initiation of transport material models 00610 Tm->initmaterialmodels(); 00611 00612 i=0; 00613 ii=0; 00614 print_init(-1, "wt"); 00615 print_step(lcid,ii,Mp->time,f); 00616 print_flush(); 00617 print_initt(-1, "wt"); 00618 print_stept(0,i,Tp->time,NULL); 00619 print_flusht(); 00620 00621 00622 // ************************** 00623 // main iteration loop *** 00624 // ************************** 00625 00626 do{ 00627 00628 // ********************************************************* 00629 // transport part 00630 // ********************************************************* 00631 00632 if (Mesprc==1) fprintf (stdout,"\n\n -----------------------------------------------------------------------------"); 00633 if (Mesprc==1) fprintf (stdout,"\n metr time %e, trfel time %e, mefel time %e",Cp->time,Tp->time,Mp->time); 00634 if (Mesprc==1) fprintf (stdout,"\n ---------------------------------------------------------------------------\n"); 00635 00636 if (Mesprc==1) fprintf (stdout,"\n\n provadi se TRFEL, trfel time %e",Tp->time); 00637 00638 Tm->updateipval (); 00639 i++; 00640 00641 // backup of attained nodal values and their time derivatives 00642 for (j=0;j<n;j++){ 00643 lhsb[j]=lhs[j]; 00644 tdlhsb[j]=tdlhs[j]; 00645 } 00646 00647 // capacity matrix 00648 capacity_matrix (lcid); 00649 00650 fprintf (Outt,"\n matice kapacity"); 00651 Cmat->printmat (Outt); 00652 00653 // conductivity matrix 00654 conductivity_matrix (lcid); 00655 00656 fprintf (Outt,"\n matice vodivosti"); 00657 Kmat->printmat (Outt); 00658 00659 //fprintf (Outt,"\n nejvetsi vl. cislo %le",Kmat->power_method (diag,1000,1.0e-6)); 00660 //fprintf (Outt,"\n nejmensi vl. cislo %le",Kmat->inverse_iteration (diag,1000,1.0e-6)); 00661 //conductivity_matrix (lcid); 00662 00663 00664 // predictor d 00665 for (j=0;j<n;j++){ 00666 d[j] = lhs[j]+dt*(1.0-alpha)*tdlhs[j]; 00667 } 00668 00669 // C.d 00670 Cmat->gmxv (d,p); 00671 00672 // matrix of the system of equations 00673 // C + alpha.dt.K 00674 Kmat->scalgm (dt*alpha); 00675 Kmat->addgm (1.0,*Cmat); 00676 00677 00678 // prescribed fluxes f_{n+1} 00679 trfel_right_hand_side (lcid,f,n); 00680 00681 // C.d + alpha dt f_n+1 00682 norrhs=0.0; 00683 for (j=0;j<n;j++){ 00684 rhs[j] = f[j]*alpha*dt + p[j]; 00685 // fprintf (stdout,"\n rhs %15.10le k je %d",rhs[j],j); 00686 //norrhs+=f[j]*f[j]; 00687 //norrhs+=rhs[j]*rhs[j]; 00688 //z[j]=rhs[j]; 00689 } 00690 //norrhs=sqrt(norrhs); 00691 00692 //Kmat->printmat (Outt); 00693 00694 // solution of the system of equations 00695 // nodal values d_{n+1} are solved 00696 00697 00698 Kmat->diag_scale (diag); 00699 for (k=0;k<n;k++){ 00700 rhs[k]*=diag[k]; 00701 } 00702 00703 Kmat->printmat (Outt); 00704 00705 Psolt->par_linear_solver (Gtt,Kmat,lhs,rhs,Outt,Mespr); 00706 00707 for (k=0;k<n;k++){ 00708 lhs[k]*=diag[k]; 00709 } 00710 00711 00712 // computation of time derivatives of the nodal values 00713 for (j=0;j<n;j++){ 00714 tdlhs[j]=(lhs[j]-d[j])/dt/alpha; 00715 } 00716 00717 // physically corrected solution 00718 solution_correction (); 00719 // approximation of nodal values into ontegration points 00720 approximation (); 00721 00722 // nulleqother (); 00723 actual_previous_nodval (); 00724 00725 00726 00727 // iteration for equilibrium 00728 for (j=0;j<ini;j++){ 00729 00730 // capacity matrix 00731 capacity_matrix (lcid); 00732 00733 // conductivity matrix 00734 conductivity_matrix (lcid); 00735 00736 00737 // matrix of the system of equations 00738 // C + alpha.dt.K 00739 Kmat->scalgm (dt*alpha); 00740 Kmat->addgm (1.0,*Cmat); 00741 00742 // (C + alpha.dt.K).d 00743 Kmat->gmxv (lhs,v); 00744 00745 // C.d 00746 Cmat->gmxv (d,p); 00747 00748 // computation of residuals 00749 //norrhs=0.0; 00750 //norres=0.0; 00751 for (k=0;k<n;k++){ 00752 rhs[k]=f[k]*alpha*dt+p[k]-v[k]; 00753 // fprintf (stdout,"\n rhs %15.10le k je %d",rhs[k],k); 00754 //norres+=rhs[k]*rhs[k]; 00755 z[k]=f[k]*alpha*dt+p[k]; 00756 //norrhs+=z[k]*z[k]; 00757 //fprintf (stdout,"\n rhs %15.10le z %15.10le v %15.10le",rhs[k],z[k],v[k]); 00758 } 00759 //norres=sqrt(norres); 00760 //norrhs=sqrt(norrhs); 00761 00762 //fprintf (stdout,"\n KONTROLA norres %15.10le",norres); 00763 00764 stop = norm_computation2 (rhs,z,err,2,1); 00765 00766 //fprintf (stdout,"\n pauza \n"); 00767 00768 if (stop==1){ 00769 break; 00770 } 00771 00772 if (Mesprt != 0) fprintf (stdout,"\n iteration number %ld",j); 00773 00774 00775 for (k=0;k<n;k++){ 00776 rhs[k]=f[k]*alpha*dt+p[k]; 00777 //z[k]=lhs[k]; 00778 } 00779 00780 Kmat->diag_scale (diag); 00781 for (k=0;k<n;k++){ 00782 rhs[k]*=diag[k]; 00783 } 00784 00785 00786 Psolt->par_linear_solver (Gtt,Kmat,lhs,rhs,Outt,Mespr); 00787 00788 for (k=0;k<n;k++){ 00789 lhs[k]*=diag[k]; 00790 } 00791 00792 // computation of time derivatives of the nodal values 00793 for (k=0;k<n;k++){ 00794 tdlhs[k]=(lhs[k]-d[k])/dt/alpha; 00795 } 00796 00797 // physically corrected solution 00798 solution_correction (); 00799 // approximation of nodal values into ontegration points 00800 approximation (); 00801 00802 00803 // nulleqother (); 00804 actual_previous_nodval (); 00805 00806 } 00807 00808 00809 // ********************************************************* 00810 // end of transport part and beginning of mechanical part 00811 // ********************************************************* 00812 00813 if (Mp->time <= Tp->time){ 00814 00815 // approximation of temperatures from TRFEL nodes to MEFEL integration points 00816 approximation_temper (); 00817 // approximation of humidities from TRFEL nodes to MEFEL integration points 00818 approximation_humid (); 00819 00820 // update of values stored at integration points 00821 Mm->updateipval(); 00822 00823 // assembling of the right hand side (prescribed load) 00824 mefel_right_hand_side (lcid,f); 00825 00826 // computation of sums of force components in particular directions 00827 // it is used in special creep and consolidation models 00828 Mb->comp_sum (f); 00829 00830 // computation of unbalanced forces which will be on the right hand side 00831 incr_internal_forces (lcid,fi); 00832 00833 for (j=0;j<nm;j++){ 00834 fb[j]=f[j]-fp[j]+fi[j]; 00835 fp[j]=f[j]; 00836 } 00837 00838 // assembling of tangent stiffness matrix 00839 stiffness_matrix (lcid); 00840 00841 // solution of K(r).v=F 00842 Psolm->par_linear_solver (Gtm,Smat,dr,fb,Out,Mespr); 00843 00844 // new displacement vector 00845 for (j=0;j<nm;j++){ 00846 r[j]+=dr[j]; 00847 } 00848 00849 // computation of internal forces 00850 internal_forces (lcid,fi); 00851 00852 // temperature contributions 00853 nullv (fb,nm); 00854 Mb->dlc[lcid].tempercontrib (lcid,fb,nm,Cp->time); 00855 00856 norfa=ss(f,f,nm); 00857 for (j=0;j<nm;j++){ 00858 fb[j] = f[j] - fi[j] - fb[j]; 00859 //fb[j] = f[j] - fi[j]; 00860 } 00861 00862 // norm of vector of unbalanced forces 00863 norf=ss(fb,fb,nm); 00864 if (norfa > 0.0) 00865 norf /= norfa; 00866 00867 00868 if (norf<err){ 00869 00870 00871 // computation of values required for output 00872 compute_req_val (lcid); 00873 // time increment 00874 Mp->time = Mp->timecon.newtime (); 00875 Mp->timecon.take_values (Mp->timecon); 00876 00877 if (Mp->timecon.isitimptime ()==1){ 00878 mtsolver_save (aux,r,fp,ii,Mp->time,nm,0,0); 00879 } 00880 00881 print_step(lcid, ii, Mp->time, f); 00882 print_flush(); 00883 ii++; 00884 } 00885 00886 // iteration of unbalanced forces caused by time independent models 00887 // internal iteration loop 00888 else{ 00889 for (j=0;j<ini;j++) 00890 { 00891 00892 if(Cp->fcsolv == fullnewtonc){ 00893 //cleaning of stiffness matrix 00894 if (Smat != NULL){ 00895 delete Smat; 00896 Smat=NULL; 00897 } 00898 //assembling of tangent stiffness matrix for fullnewton 00899 stiffness_matrix (lcid); 00900 } 00901 00902 // solution of K(r).v=F 00903 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out); 00904 00905 for (k=0;k<nm;k++) 00906 r[k]+=dr[k]; 00907 00908 // computation of internal forces 00909 internal_forces (lcid,fi); 00910 00911 // vector of unbalanced forces 00912 for (k=0;k<nm;k++) 00913 fb[k]=f[k]-fi[k]; 00914 00915 // norm of vector of unbalanced forces 00916 norf=ss(fb,fb,nm); 00917 if (norfa > 0.0) 00918 norf /= norfa; 00919 00920 if (Mespr==1) fprintf (stdout,"\n Time increment %lf inner loop j %ld norf=%e",Mp->time,j,norf); 00921 if (norf<err) break; 00922 } 00923 if (j==ini || norf>err) 00924 { 00925 if (Mespr==1) fprintf (stdout,"\n Iteration process does not converge req. err = %e, reached err = %e", err,norf); 00926 00927 // equilibrium state was not reached 00928 compute_req_val (lcid); 00929 // time increment 00930 Mp->time = Mp->timecon.newtime (); 00931 Mp->timecon.take_values (Mp->timecon); 00932 00933 if (Mp->timecon.isitimptime ()==1) 00934 mtsolver_save (aux,r,fp,ii,Mp->time,nm,0,0); 00935 print_step(lcid, ii, Mp->time, f); 00936 print_flush(); 00937 ii++; 00938 } 00939 else 00940 { 00941 // equilibrium state was reached 00942 compute_req_val (lcid); 00943 // time increment 00944 Mp->time = Mp->timecon.newtime (); 00945 Mp->timecon.take_values (Mp->timecon); 00946 00947 if (Mp->timecon.isitimptime ()==1) 00948 mtsolver_save (aux,r,fp,ii,Mp->time,nm,0,0); 00949 print_step(lcid, ii, Mp->time, f); 00950 print_flush(); 00951 ii++; 00952 } 00953 } 00954 } 00955 00956 00957 // new time and new time increment 00958 Cp->time = Cp->timecon.newtime (); 00959 00960 00961 //Mp->time = Cp->time/86400.0; 00962 //Mp->timecon.take_values (Cp->timecon); 00963 //Mp->timecon.seconds_days (); 00964 00965 00966 Mp->time = Cp->time; 00967 Tp->time = Cp->time; 00968 00969 00970 dt = Cp->timecon.actualbacktimeincr (); 00971 i++; 00972 00973 00974 }while(Cp->time<end_time); 00975 00976 print_close (); 00977 print_closet (); 00978 00979 delete [] fi; 00980 delete [] fb; 00981 delete [] fp; 00982 delete [] dr; 00983 00984 delete [] p; 00985 delete [] d; 00986 00987 delete [] fit; 00988 delete [] fbt; 00989 delete [] lhst_last; 00990 */ 00991 } 00992