source: issm/oecreview/Archive/14312-15392/ISSM-15351-15352.diff@ 15393

Last change on this file since 15393 was 15393, checked in by Mathieu Morlighem, 12 years ago

NEW: adding Archive/14312-15392 for oecreview
File size: 8.6 KB

  • TabularUnified ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

     
    88#include "../modules/modules.h"
    99
    1010void solutionsequence_hydro_nonlinear(FemModel* femmodel){
    11 
    1211        /*solution : */
    1312        Vector<IssmDouble>* ug_sed=NULL;
    1413        Vector<IssmDouble>* ug_epl=NULL;
    1514        Vector<IssmDouble>* uf=NULL;
    16         Vector<IssmDouble>* old_uf=NULL;
    17         Vector<IssmDouble>* aged_ug_sed=NULL;
    18         Vector<IssmDouble>* aged_ug_epl=NULL;
     15        Vector<IssmDouble>* uf_int_iter=NULL;
     16        Vector<IssmDouble>* ug_sed_main_iter=NULL;
     17        Vector<IssmDouble>* ug_epl_main_iter=NULL;
    1918        Vector<IssmDouble>* ys=NULL;
    2019        Vector<IssmDouble>* dug=NULL;
    21 
     20        Vector<IssmDouble>* old_ug=NULL;
     21       
    2222        Matrix<IssmDouble>* Kff=NULL;
    2323        Matrix<IssmDouble>* Kfs=NULL;
    2424        Vector<IssmDouble>* pf=NULL;
    … …  
    3333        IssmDouble sediment_kmax,time;
    3434        IssmDouble eps_hyd;
    3535        IssmDouble ndu_sed,nu_sed;
    36         IssmDouble ndu_epl,nu_epl;     
     36        IssmDouble ndu_epl,nu_epl;
    3737
    3838        /*Recover parameters: */
    3939        femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME
    4040        femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    4141        femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum);
    4242        femmodel->parameters->FindParam(&time,TimeEnum);
    43         hydro_maxiter=100;
     43        hydro_maxiter=150;
    4444        hydrocount=1;
    4545
    4646        /*Iteration on the two layers*/
    … …  
    5656                sedcount=1;
    5757                eplcount=1;
    5858                //save pointer to old velocity
    59                 delete aged_ug_sed;
    60                 aged_ug_sed=ug_sed;
     59                delete ug_sed_main_iter;
     60                ug_sed_main_iter=ug_sed;
    6161                if(isefficientlayer){
    62                         delete aged_ug_epl;
    63                         aged_ug_epl=ug_epl;
     62                        delete ug_epl_main_iter;
     63                        ug_epl_main_iter=ug_epl;
    6464                }
    6565
    6666                femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
    … …  
    7777                        CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
    7878                        Reduceloadx(pf,Kfs,ys); delete Kfs;
    7979                        delete uf;
    80                         Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
    81                         delete old_uf; old_uf=uf->Duplicate();
    82                         if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting aged_ug_sed*/
     80                        Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
     81                        delete uf_int_iter; uf_int_iter=uf->Duplicate();
     82                        if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting ug_sed_main_iter*/
    8383                        delete Kff; delete pf; delete df;
    8484
    8585                        Mergesolutionfromftogx(&ug_sed,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
    … …  
    119119                                CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
    120120                                Reduceloadx(pf,Kfs,ys); delete Kfs;
    121121                                delete uf;
    122                                 Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
    123                                 delete old_uf; old_uf=uf->Duplicate();
     122                                Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
     123                                delete uf_int_iter; uf_int_iter=uf->Duplicate();
    124124                                if(eplcount>1) delete ug_epl;
    125125                                delete Kff;delete pf;
    126126                                delete df;
    127127                                Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
    128128                                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
    129                                 //femmodel->HydrologyEPLupdateDomainx();
    130129                                ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     130                                femmodel->HydrologyEPLupdateDomainx();                 
    131131                               
    132132                                if (!eplconverged){
    133133                                        if(VerboseConvergence()) _printf0_("   # EPL unstable constraints = " << num_unstable_constraints << "\n");
    … …  
    137137                                        }
    138138                                }
    139139                                eplcount++;
    140                                
     140
    141141                                if(eplconverged){
    142142                                        InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum);
    143143                                        InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum);
    144144                                        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
    145                                         //femmodel->HydrologyEPLupdateDomainx();
    146145                                        break;
    147146                                }
    148147                        }
    … …  
    151150                if(!hydroconverged){
    152151                        //compute norm(du)/norm(u)
    153152                        dug=ug_sed->Duplicate(); _assert_(dug);
    154                         aged_ug_sed->Copy(dug);
     153                        ug_sed_main_iter->Copy(dug);   
    155154                        dug->AYPX(ug_sed,-1.0);
    156                         ndu_sed=dug->Norm(NORM_TWO); nu_sed=aged_ug_sed->Norm(NORM_TWO);
     155                        ndu_sed=dug->Norm(NORM_TWO); nu_sed=ug_sed_main_iter->Norm(NORM_TWO);
    157156                        if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
    158157                        if (!xIsNan<IssmDouble>(eps_hyd)){
    159158                                if (!isefficientlayer){
    160159                                        if ((ndu_sed/nu_sed)<eps_hyd){
    161                                                 if(VerboseConvergence()) _printf0_(setw(50) << left << "   Converged \n");
     160                                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Converged after, " << hydrocount << " iterations \n");
    162161                                                hydroconverged=true;
    163162                                        }
    164163                                        else{
    … …  
    168167                                }
    169168                                else{
    170169                                        dug=ug_epl->Duplicate();_assert_(dug);
    171                                         aged_ug_epl->Copy(dug);_assert_(aged_ug_epl);
     170                                        ug_epl_main_iter->Copy(dug);_assert_(ug_epl_main_iter);
    172171                                        dug->AYPX(ug_epl,-1.0);
    173172                                        ndu_epl=dug->Norm(NORM_TWO);
    174                                         nu_epl=aged_ug_epl->Norm(NORM_TWO);
     173                                        nu_epl=ug_epl_main_iter->Norm(NORM_TWO);
    175174
    176175                                        if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!");
    177176                                        if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
    178                                         if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<eps_hyd){
    179                                                 if (VerboseConvergence()) _printf0_(setw(50) << left << "   Converged \n");
     177                                        if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<(eps_hyd*10)){
     178                                                if (VerboseConvergence()) _printf0_(setw(50) << left << "   Converged after, " << hydrocount << " iterations \n");
    180179                                                hydroconverged=true;
    181180                                        }
    182181                                        else{
    … …  
    201200        delete ug_epl;
    202201        delete ug_sed;
    203202        delete uf;
    204         delete old_uf;
    205         delete aged_ug_sed;
    206         delete aged_ug_epl;
     203        delete uf_int_iter;
     204        delete ug_sed_main_iter;
     205        delete ug_epl_main_iter;
    207206        delete dug;
    208 
    209207}

  • TabularUnified ../trunk-jpl/src/c/classes/FemModel.cpp

     
    852852                        element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    853853                        element->CreateKMatrix(Kff_temp,NULL);
    854854                }
     855
    855856                for (i=0;i<this->loads->Size();i++){
    856857                        load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    857858                        if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
    … …  
    871872                element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    872873                element->CreateKMatrix(Kff,Kfs);
    873874        }
     875       
    874876        for (i=0;i<this->loads->Size();i++){
    875877                load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    876878                if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);

  • TabularUnified ../trunk-jpl/src/c/classes/Elements/Tria.cpp

     
    61316131        for(int ig=gauss->begin();ig<gauss->end();ig++){
    61326132
    61336133                gauss->GaussPoint(ig);
    6134 
    61356134                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    61366135
    61376136                /*Diffusivity*/
    … …  
    61566155                                                &Ke->values[0],1);
    61576156                }
    61586157        }
    6159 
    61606158        /*Clean up and return*/
    61616159        delete gauss;
    61626160        return Ke;
    … …  
    62436241
    62446242        /* Start  looping on the number of gaussian points: */
    62456243        gauss=new GaussTria(2);
    6246 
    62476244        for(int ig=gauss->begin();ig<gauss->end();ig++){
    62486245
    62496246                gauss->GaussPoint(ig);
    6250 
    62516247                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    62526248                GetNodalFunctions(basis, gauss);
    62536249
    … …  
    63106306        for(int ig=gauss->begin();ig<gauss->end();ig++){
    63116307
    63126308                gauss->GaussPoint(ig);
    6313 
    63146309                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    63156310                GetNodalFunctions(basis, gauss);
    63166311
    … …  
    66906685
    66916686                /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
    66926687                this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
    6693                 if(eplhead[i]>=h_max && this->AnyActive()){                     
     6688                if(eplhead[i]>=h_max && this->AnyActive()){
    66946689                        for(j=0;j<numdof;j++){
    66956690                                if(old_active[j]>0.){
    66966691                                        vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
Note: See TracBrowser for help on using the repository browser.