Changeset 16720


Ignore:
Timestamp:
11/12/13 15:54:29 (11 years ago)
Author:
seroussi
Message:

NEW: adding mass transport in analyses

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r16684 r16720  
    216216}/*}}}*/
    217217void MasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    218         _error_("not implemented yet");
    219 }/*}}}*/
     218
     219        int        i,hydroadjustment;
     220        int*       doflist=NULL;
     221        IssmDouble rho_ice,rho_water,minthickness;
     222        //Element* basalelement=element->SpawnBasalElement();
     223        Element* basalelement=element;
     224
     225        /*Fetch number of nodes for this finite element*/
     226        int numnodes = basalelement->GetNumberOfNodes();
     227
     228        /*Fetch dof list and allocate solution vector*/
     229        basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     230        IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
     231        IssmDouble* newbed       = xNew<IssmDouble>(numnodes);
     232        IssmDouble* newsurface   = xNew<IssmDouble>(numnodes);
     233        IssmDouble* oldthickness = xNew<IssmDouble>(numnodes);
     234        IssmDouble* oldbed       = xNew<IssmDouble>(numnodes);
     235        IssmDouble* oldsurface   = xNew<IssmDouble>(numnodes);
     236        IssmDouble* phi          = xNew<IssmDouble>(numnodes);
     237
     238        /*Use the dof list to index into the solution vector: */
     239        basalelement->FindParam(&minthickness,MasstransportMinThicknessEnum);
     240        for(i=0;i<numnodes;i++){
     241                newthickness[i]=solution[doflist[i]];
     242                if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
     243                /*Constrain thickness to be at least 1m*/
     244                if(newthickness[i]<minthickness) newthickness[i]=minthickness;
     245        }
     246
     247        /*Get previous bed, thickness and surface*/
     248        basalelement->GetInputListOnNodes(&oldbed[0],BedEnum);
     249        basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
     250        basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
     251        basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
     252
     253        /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
     254        basalelement->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
     255        rho_ice   = basalelement->GetMaterialParameter(MaterialsRhoIceEnum);
     256        rho_water = basalelement->GetMaterialParameter(MaterialsRhoWaterEnum);
     257
     258        for(i=0;i<numnodes;i++) {
     259                /*If shelf: hydrostatic equilibrium*/
     260                if (phi[i]>0.){
     261                        newsurface[i] = oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
     262                        newbed[i]     = oldbed[i];                 //same bed: do nothing
     263                }
     264                else{ //this is an ice shelf
     265                        if(hydroadjustment==AbsoluteEnum){
     266                                newsurface[i] = newthickness[i]*(1-rho_ice/rho_water);
     267                                newbed[i]     = newthickness[i]*(-rho_ice/rho_water);
     268                        }
     269                        else if(hydroadjustment==IncrementalEnum){
     270                                newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
     271                                newbed[i]     = oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed               = oldbed + di * dH
     272                        }
     273                        else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
     274                }
     275        }
     276
     277        /*Add input to the element: */
     278        element->AddBasalInput(ThicknessEnum,newthickness,P1Enum);
     279        element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
     280        element->AddBasalInput(BedEnum,newbed,P1Enum);
     281
     282        /*Free ressources:*/
     283        xDelete<IssmDouble>(newthickness);
     284        xDelete<IssmDouble>(newbed);
     285        xDelete<IssmDouble>(newsurface);
     286        xDelete<IssmDouble>(oldthickness);
     287        xDelete<IssmDouble>(oldbed);
     288        xDelete<IssmDouble>(oldsurface);
     289        xDelete<IssmDouble>(phi);
     290        xDelete<int>(doflist);
     291        if(basalelement!=element) delete element;
     292}/*}}}*/

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

    r16717 r16720  
    3939                virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;
    4040                virtual void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
     41                virtual void   AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
    4142                virtual void   AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
    4243                virtual void   AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;

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

    r16716 r16720  
    126126        _assert_(this->inputs);
    127127        this->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
     128}
     129/*}}}*/
     130/*FUNCTION Penta::AddBasalInput{{{*/
     131void  Penta::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){
     132
     133        _assert_(this->inputs);
     134        if(!IsOnBed()) return;
     135        else{
     136                if(interpolation_enum==P1Enum){
     137                        int        i;
     138                        IssmDouble extrudedvalues[NUMVERTICES];
     139                        Penta*     penta=NULL;
     140
     141                        for(i=1;i<NUMVERTICES2D;i++){
     142                                extrudedvalues[i]=values[i];
     143                                extrudedvalues[i+NUMVERTICES2D]=values[i];
     144                        }
     145                        this->inputs->AddInput(new PentaInput(input_enum,&extrudedvalues[0],P1Enum));
     146                        penta=this;
     147                        for(;;){
     148                                penta->inputs->AddInput(new PentaInput(input_enum,&extrudedvalues[0],P1Enum));
     149                                if (penta->IsOnSurface()) break;
     150                                penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
     151                        }
     152                }
     153                else _error_("not implemented yet");
     154        }
    128155}
    129156/*}}}*/

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

    r16715 r16720  
    184184                /*}}}*/
    185185                /*Penta specific routines:{{{*/
     186                void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
    186187                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
    187188                void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r16717 r16720  
    6868                /*}}}*/
    6969                /*Element virtual functions definitions: {{{*/
     70                void        AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
    7071                void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
    7172                void        AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};

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

    r16716 r16720  
    164164        *pd_nz=d_nz;
    165165        *po_nz=o_nz;
     166}
     167/*}}}*/
     168/*FUNCTION Tria::AddBasalInput{{{*/
     169void  Tria::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){
     170
     171        /*Call inputs method*/
     172        _assert_(this->inputs);
     173       
     174        int meshtype;
     175        parameters->FindParam(&meshtype,MeshTypeEnum);
     176        switch(meshtype){
     177                case Mesh2DhorizontalEnum:
     178                        this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
     179                        break;
     180                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     181        }
     182
    166183}
    167184/*}}}*/

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

    r16715 r16720  
    195195                /*}}}*/
    196196                /*Tria specific routines:{{{*/
     197                void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
    197198                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
    198199                void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
Note: See TracChangeset for help on using the changeset viewer.