source: issm/oecreview/Archive/14312-15392/ISSM-15222-15223.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: 17.2 KB

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

     
    1010void solutionsequence_hydro_nonlinear(FemModel* femmodel){
    1111
    1212        /*solution : */
    13         Vector<IssmDouble>* ug=NULL;
    14         Vector<IssmDouble>* uf_sed=NULL;
    15         Vector<IssmDouble>* uf_epl=NULL;
     13        Vector<IssmDouble>* ug_sed=NULL;
     14        Vector<IssmDouble>* ug_epl=NULL;
     15        Vector<IssmDouble>* uf=NULL;
    1616        Vector<IssmDouble>* old_uf=NULL;
    17         Vector<IssmDouble>* aged_uf_sed=NULL;
    18         Vector<IssmDouble>* aged_uf_epl=NULL;
     17        Vector<IssmDouble>* aged_ug_sed=NULL;
     18        Vector<IssmDouble>* aged_ug_epl=NULL;
    1919        Vector<IssmDouble>* ys=NULL;
    20         Vector<IssmDouble>* duf=NULL;
     20        Vector<IssmDouble>* dug=NULL;
    2121
    2222        Matrix<IssmDouble>* Kff=NULL;
    2323        Matrix<IssmDouble>* Kfs=NULL;
    … …  
    4949
    5050        /*Iteration on the two layers + transfer*/
    5151        femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
    52         GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    53         Reducevectorgtofx(&uf_sed, ug, femmodel->nodes,femmodel->parameters); _assert_(uf_sed);
     52        GetSolutionFromInputsx(&ug_sed, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    5453        if(isefficientlayer) {
    5554                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
    56                 GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    57                 Reducevectorgtofx(&uf_epl, ug, femmodel->nodes,femmodel->parameters);
    58         }
     55                GetSolutionFromInputsx(&ug_epl, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
     56}
    5957
    6058        for(;;){
    6159                sedcount=1;
    6260                eplcount=1;
    6361                //save pointer to old velocity
    64                 delete aged_uf_sed;aged_uf_sed=uf_sed;
     62                delete aged_ug_sed;
     63                aged_ug_sed=ug_sed;
    6564                if(isefficientlayer){
    66                         delete aged_uf_epl;
    67                         aged_uf_epl=uf_epl;
     65                        delete aged_ug_epl;
     66                        aged_ug_epl=ug_epl;
    6867                }
    6968
    7069                femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
    … …  
    7271                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum);
    7372                femmodel->UpdateConstraintsx();
    7473                femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
    75                 /*              femmodel->HydrologyTransferx();
    76                  */
     74
    7775                /*Iteration on the sediment layer*/
    7876                for(;;){
    7977                        femmodel->HydrologyTransferx();
    8078                        femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax);
    8179                        CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
    8280                        Reduceloadx(pf,Kfs,ys); delete Kfs;
    83                         if(sedcount>1)delete uf_sed;
    84                         Solverx(&uf_sed, Kff, pf,old_uf, df, femmodel->parameters);
    85                         delete old_uf; old_uf=uf_sed->Duplicate();
    86                         delete Kff; delete pf; delete ug; delete df;
     81                        delete uf;
     82                        Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
     83                        delete old_uf; old_uf=uf->Duplicate();
     84                        if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting aged_ug_sed*/
     85                        delete Kff; delete pf; delete df;
    8786
    88                         Mergesolutionfromftogx(&ug,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys;
    89                         InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
     87                        Mergesolutionfromftogx(&ug_sed,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
     88                        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_sed);
    9089                        ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    9190
    9291                        if (!sedconverged){
    … …  
    101100                        if(sedconverged){
    102101                                femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
    103102                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sedconverged,ConvergedEnum);
    104                                 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
     103                                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_sed);
    105104                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);
    106105                                break;
    107106                        }
    … …  
    121120                        femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL);
    122121                        CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
    123122                        Reduceloadx(pf,Kfs,ys); delete Kfs;
    124                         if(eplcount>1) delete uf_epl;
    125                         Solverx(&uf_epl, Kff, pf,old_uf, df, femmodel->parameters);
    126                         delete old_uf; old_uf=uf_epl->Duplicate();
    127                         delete Kff;delete pf; delete ug; delete df;
    128                         Mergesolutionfromftogx(&ug,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
    129                         InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
    130 
     123                        delete uf;
     124                        Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
     125                        delete old_uf; old_uf=uf->Duplicate();                 
     126                        if(eplcount>1) delete ug_epl;
     127                        delete Kff;delete pf;
     128                        delete df;
     129                        Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
     130                        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
    131131                        ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    132132
    133133                        if (!eplconverged){
    … …  
    142142                        if(eplconverged){
    143143                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum);
    144144                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum);
    145                                 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
     145                                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
    146146                                break;
    147147                        }
    148148                }
    … …  
    150150                /*System convergence check*/
    151151                if(!hydroconverged){
    152152                        //compute norm(du)/norm(u)
    153                         _assert_(aged_uf_sed); _assert_(uf_sed);
    154                         duf=uf_sed->Duplicate(); _assert_(duf);
    155                         aged_uf_sed->Copy(duf);
    156                         duf->AYPX(uf_sed,-1.0);
    157                         ndu_sed=duf->Norm(NORM_TWO); nu_sed=aged_uf_sed->Norm(NORM_TWO);
     153                        dug=ug_sed->Duplicate(); _assert_(dug);
     154                        aged_ug_sed->Copy(dug);
     155                        dug->AYPX(ug_sed,-1.0);
     156                        ndu_sed=dug->Norm(NORM_TWO); nu_sed=aged_ug_sed->Norm(NORM_TWO);
    158157                        if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
    159158                        if (!xIsNan<IssmDouble>(eps_hyd)){
    160159                                if (!isefficientlayer){
    … …  
    168167                                        }
    169168                                }
    170169                                else{
    171                                         duf=aged_uf_epl->Duplicate(); aged_uf_epl->Copy(duf); duf->AYPX(uf_epl,-1.0);
    172                                         ndu_epl=duf->Norm(NORM_TWO); nu_epl=aged_uf_epl->Norm(NORM_TWO);
     170                                        dug=ug_epl->Duplicate();_assert_(dug);
     171                                        aged_ug_epl->Copy(dug);_assert_(aged_ug_epl);
     172                                        dug->AYPX(ug_epl,-1.0);
     173                                        ndu_epl=dug->Norm(NORM_TWO);
     174                                        nu_epl=aged_ug_epl->Norm(NORM_TWO);
     175
    173176                                        if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!");
    174177                                        if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
    175178                                        if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<eps_hyd){
    … …  
    192195                if(hydroconverged)break;
    193196        }
    194197
    195         InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
     198        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_sed);
     199        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_epl);
    196200       
    197201        /*Free ressources: */
    198         delete ug;
    199         delete uf_sed;
    200         delete uf_epl;
     202        delete ug_epl;
     203        delete ug_sed;
     204        delete uf;
    201205        delete old_uf;
    202         delete aged_uf_sed;
    203         delete aged_uf_epl;
    204         delete duf;
     206        delete aged_ug_sed;
     207        delete aged_ug_epl;
     208        delete dug;
    205209
    206210}

  • ../trunk-jpl/src/c/classes/DofIndexing.cpp

     
    159159        else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
    160160}
    161161/*}}}*/
    162 /*FUNCTION DofIndexing::Deactivate{{{*/
    163 void DofIndexing::Deactivate(void){
    164         this->active = false;
    165162
    166         /*Constrain to 0. at this point*/
    167         for(int i=0;i<this->gsize;i++){
    168                 this->f_set[i]    = false;
    169                 this->s_set[i]    = true;
    170                 this->svalues[i]  = 0.;
    171         }
    172         return;
    173 }
    174 /*}}}*/
    175 
    176163/*Some of the Object functionality: */
    177164/*FUNCTION DofIndexing::Echo{{{*/
    178165void DofIndexing::Echo(void){
    … …  
    234221        _printf_("\n");
    235222}               
    236223/*}}}*/
     224/*FUNCTION DofIndexing::Deactivate{{{*/
     225void DofIndexing::Deactivate(void){
     226        this->active = false;
     227
     228        /*Constrain to 0. at this point*/
     229        for(int i=0;i<this->gsize;i++){
     230                this->f_set[i]    = false;
     231                this->s_set[i]    = true;
     232                this->svalues[i]  = 0.;
     233        }
     234        return;
     235}
     236/*}}}*/
     237/*FUNCTION DofIndexing::Activate{{{*/
     238void DofIndexing::Activate(void){
     239        this->active = true;
     240
     241        /*Constrain to 0. at this point*/
     242        for(int i=0;i<this->gsize;i++){
     243                this->f_set[i]    = true;
     244                this->s_set[i]    = false;
     245                this->svalues[i]  = 0.;
     246        }
     247        return;
     248}
     249/*}}}*/

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

     
    667667/*}}}*/
    668668void FemModel::UpdateConstraintsx(void){ /*{{{*/
    669669
     670        Element   *element = NULL;
    670671        IssmDouble time;
    671         int    analysis_type;
     672        int        analysis_type;
    672673
    673674        /*retrieve parameters: */
    674675        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    677678        /*start module: */
    678679        if(VerboseModule()) _printf0_("   Updating constraints for time: " << time << "\n");
    679680
    680         /*First, update dof constraints in nodes, using constraints: */
     681        /*First, Nodes might be activated/deactivated by element*/
     682        for(int i=0;i<elements->Size();i++){
     683                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     684                element->UpdateConstraints();
     685        }
     686
     687        /*Second, constraints might be time dependent: */
    681688        SpcNodesx(nodes,constraints,parameters,analysis_type);
    682689
    683690        /*Now, update degrees of freedoms: */

  • ../trunk-jpl/src/c/classes/Node.h

     
    8686                void   GetDofList(int* poutdoflist,int approximation_enum,int setenum);
    8787                void   GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum);
    8888                void   FreezeDof(int dof);
     89                void   Activate(void);
    8990                void   Deactivate(void);
    9091                int    IsFloating();
    9192                int    IsGrounded();

  • ../trunk-jpl/src/c/classes/Elements/Element.h

     
    134134                #ifdef _HAVE_HYDROLOGY_
    135135                virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
    136136                virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
     137                virtual void UpdateConstraints(void)=0;
    137138                #endif
    138139};
    139140#endif

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

     
    28192819        this->parameters=NULL;
    28202820}
    28212821/*}}}*/
     2822/*FUNCTION Tria::UpdateConstraints{{{*/
     2823void  Tria::UpdateConstraints(void){
     2824
     2825        int analysis_type;
     2826        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     2827
     2828        /*Skip if water element*/
     2829        if(IsOnWater()) return;
     2830
     2831        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     2832        switch(analysis_type){
     2833                #ifdef _HAVE_HYDROLOGY_
     2834                case HydrologyDCEfficientAnalysisEnum:
     2835                        this->HydrologyUpdateEplDomain();
     2836                        break;
     2837                #endif
     2838        }
     2839
     2840}
     2841/*}}}*/
    28222842/*FUNCTION Tria::UpdatePotentialUngrounding{{{*/
    28232843int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
    28242844
    … …  
    64066426        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    64076427
    64086428        /*Get the flag to know if the efficient layer is present*/
    6409                 this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     6429        this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    64106430
    64116431        /*Use the dof list to index into the solution vector: */
    64126432        for(int i=0;i<numdof;i++){
    … …  
    65906610}
    65916611
    65926612/*}}}*/
     6613/*FUNCTION Tria::HydrologyUpdateEplDomain{{{*/
     6614void  Tria::HydrologyUpdateEplDomain(void){
     6615
     6616        /*Intermediaries*/
     6617        const int  numdof = NDOF1 *NUMVERTICES;
     6618        bool       isefficientlayer;
     6619        IssmDouble activeEpl[numdof];
     6620
     6621        /*Get parameter and Inout list*/
     6622        this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     6623        GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
     6624
     6625        if(isefficientlayer){
     6626                for(int i=0;i<numdof;i++){
     6627                        if(activeEpl[i]==1.0)
     6628                         nodes[i]->Activate();
     6629                        else if (activeEpl[i]==0.0)
     6630                         nodes[i]->Deactivate();
     6631                        else
     6632                         _error_("activation value "<< activeEpl[i] <<" not supported");
     6633                }
     6634        }
     6635}
     6636/*}}}*/
    65936637#endif
    65946638
    65956639#ifdef _HAVE_DAKOTA_

  • ../trunk-jpl/src/c/classes/Elements/Tria.h

     
    249249                void      GetSolutionFromInputsHydrologyDCInefficient(Vector<IssmDouble>* solution);
    250250                void      GetSolutionFromInputsHydrologyDCEfficient(Vector<IssmDouble>* solution);
    251251                void    CreateHydrologyWaterVelocityInput(void);
     252                void    HydrologyUpdateEplDomain(void);
     253                void    UpdateConstraints(void);
    252254                void      InputUpdateFromSolutionHydrology(IssmDouble* solution);
    253255                void      InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
    254256                void    InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);

  • ../trunk-jpl/src/c/classes/Elements/Penta.cpp

     
    31883188        return nflipped;
    31893189}
    31903190/*}}}*/
     3191/*FUNCTION Penta::UpdateConstraints{{{*/
     3192void Penta::UpdateConstraints(void){
     3193
     3194        /*Do nothing for now*/
     3195
     3196}
     3197/*}}}*/
    31913198/*FUNCTION Penta::ViscousHeatingCreateInput {{{*/
    31923199void Penta::ViscousHeatingCreateInput(void){
    31933200
    … …  
    92309237}
    92319238/*}}}*/
    92329239#endif
    9233 
    92349240#ifdef _HAVE_BALANCED_
    92359241/*FUNCTION Penta::CreateKMatrixBalancethickness {{{*/
    92369242ElementMatrix* Penta::CreateKMatrixBalancethickness(void){

  • ../trunk-jpl/src/c/classes/Elements/Penta.h

     
    303303                #endif
    304304
    305305                #ifdef _HAVE_HYDROLOGY_
    306                 void    CreateHydrologyWaterVelocityInput(void);
     306                void    UpdateConstraints(void);
    307307                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    308308                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    309309                #endif

  • ../trunk-jpl/src/c/classes/DofIndexing.h

     
    4848                /*}}}*/
    4949                /*DofIndexing management: {{{*/
    5050                DofIndexing* Spawn(int* indices, int numindices);
     51                void Activate(void);
    5152                void Deactivate(void);
    5253                /*}}}*/
    5354

  • ../trunk-jpl/src/c/classes/Node.cpp

     
    498498
    499499}
    500500/*}}}*/
     501/*FUNCTION Node::Activate{{{*/
     502void  Node::Activate(void){
     503
     504        indexing.Activate();
     505
     506}
     507/*}}}*/
    501508/*FUNCTION Node::GetApproximation {{{*/
    502509int   Node::GetApproximation(){
    503510
Note: See TracBrowser for help on using the repository browser.