Changeset 21311

Timestamp:

10/25/16 11:44:04 (8 years ago)

Author:

Eric.Larour

Message:

CHG: new time dependent GetVectorFromInputsx routine. New EarthMassBalance module in sealevelrise_core
that relies on this routine.

Location:

issm/branches/trunk-larour-NatGeoScience2016/src/c

Files:

• classes/Elements/Element.cpp (modified) (2 diffs)
• classes/Elements/Element.h (modified) (2 diffs)
• cores/cores.h (modified) (1 diff)
• cores/sealevelrise_core.cpp (modified) (3 diffs)
• modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (modified) (1 diff)
• modules/GetVectorFromInputsx/GetVectorFromInputsx.h (modified) (1 diff)

issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp

@@ -963 +963 @@
 }
 /*}}}*/
+void       Element::GetInputListOnVerticesAtTime(IssmDouble* pvalue, int enumtype, IssmDouble time){/*{{{*/
+
+        /*Recover input*/
+        Input* input=this->GetInput(enumtype);
+        if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
+
+        /*Fetch number vertices for this element*/
+        int numvertices = this->GetNumberOfVertices();
+
+        /*Checks in debugging mode*/
+        _assert_(pvalue);
+
+        /* Start looping on the number of vertices: */
+        Gauss*gauss=this->NewGauss();
+        for(int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+                input->GetInputValue(&pvalue[iv],gauss,time);
+        }
+
+        /*clean-up*/
+        delete gauss;
+}
+/*}}}*/
 void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
 
@@ -1172 +1195 @@
                 this->GetInputListOnNodes(values,input_enum);
                 vector->SetValues(numnodes,doflist,values,INS_VAL);
+                break;
+        default:
+                _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
+        }
+        
+        /*Clean up*/
+        xDelete<int>(doflist);
+        xDelete<IssmDouble>(values);
+
+}
+/*}}}*/
+void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type, IssmDouble time){/*{{{*/
+
+        /*Fetch number vertices for this element and allocate arrays*/
+        int         numvertices = this->GetNumberOfVertices();
+        int         numnodes    = this->GetNumberOfNodes();
+        int*        doflist     = NULL;
+        IssmDouble* values      = NULL;
+
+        switch(type){
+        case VertexPIdEnum:
+                doflist = xNew<int>(numvertices);
+                values = xNew<IssmDouble>(numvertices);
+                /*Fill in values*/
+                this->GetVertexPidList(doflist);
+                this->GetInputListOnVerticesAtTime(values,input_enum,time);
+                vector->SetValues(numvertices,doflist,values,INS_VAL);
                 break;
         default:

issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.h

@@ -89 +89 @@
                 void               GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype);
                 void               GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
+                void               GetInputListOnVerticesAtTime(IssmDouble* pvalue,int enumtype,IssmDouble time);
                 void               GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
                 void               GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug);
@@ -100 +101 @@
                 void               GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
                 void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum, int type);
+                void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum, int type,IssmDouble time);
                 void                 GetVertexPidList(int* pidlist);
                 void               GetVerticesConnectivityList(int* connectivitylist);

issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h

@@ -64 +64 @@
 void TransferForcing(FemModel* femmodel,int forcingenum);
 void TransferSealevel(FemModel* femmodel,int forcingenum);
+void EarthMassTransport(FemModel* femmodel);
 
 //solution configuration

issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp

@@ -40 +40 @@
         femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
         femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
-
+        
         /*first, recover lat,long and radius vectors from vertices: */
         VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
@@ -71 +71 @@
         /*set configuration: */
         if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
+
+        /*figure out whether deltathickness was computed on the earth (mass trasnport module should 
+         * have taken care of it for ice caps:  */
+        if(iscoupler){
+                int ismasstransport,modelid,earthid;
+                femmodel->parameters->FindParam(&modelid,ModelIdEnum);
+                femmodel->parameters->FindParam(&earthid,EarthIdEnum);
+                femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
+                if((modelid==EarthIdEnum) & !ismasstransport) EarthMassTransport(femmodel);
+        }
 
         /*transfer deltathickness forcing from ice caps to earth model: */
@@ -359 +369 @@
 
 } /*}}}*/
+void EarthMassTransport(FemModel* femmodel){ /*{{{*/
+
+        IssmDouble time,dt;
+        Vector<IssmDouble> *oldthickness    = NULL;
+        Vector<IssmDouble> *newthickness    = NULL;
+        Vector<IssmDouble> *deltathickness    = NULL;
+        int nv;
+
+        /*No mass transport module was called, so we are just going to retrieve the geometry thickness 
+         * at this time step, at prior time step, and plug the difference as deltathickness: */
+        femmodel->parameters->FindParam(&time,TimeEnum);
+        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+        nv=femmodel->vertices->NumberOfVertices();
+
+        GetVectorFromInputsx(&newthickness,femmodel,ThicknessEnum,VertexSIdEnum);
+        GetVectorFromInputsx(&oldthickness,femmodel,ThicknessEnum,VertexSIdEnum,time-dt);
+
+        /*compute deltathickness: */
+        deltathickness = new Vector<IssmDouble>(nv); 
+        newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); 
+                
+        /*plug into elements:*/
+        InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum);
+
+        /*free ressources:*/
+        delete oldthickness;
+        delete newthickness;
+        delete deltathickness;
+
+} /*}}}*/

issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp

@@ -33 +33 @@
         *pvector=vector;
 }
+void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time){
+
+        int i;
+        Vector<IssmDouble>* vector=NULL;
+
+        switch(type){
+        case VertexPIdEnum: case VertexSIdEnum:
+                vector=new Vector<IssmDouble>(femmodel->vertices->NumberOfVertices());
+                break;
+        case NodesEnum:case NodeSIdEnum:
+                vector=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
+                break;
+        default:
+                        _error_("vector type: " << EnumToStringx(type) << " not supported yet!");
+        }
+        /*Look up in elements*/
+        for(i=0;i<femmodel->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                element->GetVectorFromInputs(vector,name,type,time);
+        }
+
+        vector->Assemble();
+
+        /*Assign output pointers:*/
+        *pvector=vector;
+}
 
 void GetVectorFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name, int type){

issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h

@@ -9 +9 @@
 /* local prototypes: */
 void    GetVectorFromInputsx( Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type);
+void    GetVectorFromInputsx( Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time);
 void    GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type);