Changeset 16903

Timestamp:

11/23/13 08:47:33 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: done with hydrology Shreve

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyShreveAnalysis.cpp (modified) (3 diffs)
• analyses/HydrologyShreveAnalysis.h (modified) (1 diff)
• analyses/MasstransportAnalysis.cpp (modified) (3 diffs)
• classes/Materials/Matpar.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp ¶

@@ -83 +83 @@
 /*Finite Element Analysis*/
 ElementMatrix* HydrologyShreveAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        /*Intermediaries */
+        IssmDouble diffusivity;
+        IssmDouble Jdet,D_scalar,dt,h;
+        IssmDouble vx,vy,vel,dvxdx,dvydy;
+        IssmDouble dvx[2],dvy[2];
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke     = element->NewElementMatrix();
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+        IssmDouble     D[2][2]={0.};
+
+        /*Create water velocity vx and vy from current inputs*/
+        CreateHydrologyWaterVelocityInput(element);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        Input* vx_input=element->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
+        Input* vy_input=element->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
+        h = element->CharacteristicLength();
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+
+                D_scalar=gauss->weight*Jdet;
+
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D_scalar,1,1,0,
+                                        basis,1,numnodes,0,
+                                        Ke->values,1);
+
+                GetB(B,element,xyz_list,gauss);
+                GetBprime(Bprime,element,xyz_list,gauss);
+
+                dvxdx=dvx[0];
+                dvydy=dvy[1];
+                D_scalar=dt*gauss->weight*Jdet;
+
+                D[0][0]=D_scalar*dvxdx;
+                D[1][1]=D_scalar*dvydy;
+                TripleMultiply(B,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        B,2,numnodes,0,
+                                        &Ke->values[0],1);
+
+                D[0][0]=D_scalar*vx;
+                D[1][1]=D_scalar*vy;
+                TripleMultiply(B,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        Bprime,2,numnodes,0,
+                                        &Ke->values[0],1);
+
+                /*Artificial diffusivity*/
+                vel=sqrt(vx*vx+vy*vy);
+                D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
+                D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
+                D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
+                D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
+                TripleMultiply(Bprime,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        Bprime,2,numnodes,0,
+                                        &Ke->values[0],1);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        delete gauss;
+        return Ke;
 }/*}}}*/
 ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
@@ -134 +221 @@
         return pe;
 }/*}}}*/
+void HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+         * For node i, Bi can be expressed in the actual coordinate system
+         * by: 
+         *       Bi=[ N ]
+         *          [ N ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions*/
+        IssmDouble* basis=xNew<IssmDouble>(numnodes);
+        element->NodalFunctions(basis,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B[numnodes*0+i] = basis[i];
+                B[numnodes*1+i] = basis[i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(basis);
+}/*}}}*/
+void HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+         * For node i, Bi' can be expressed in the actual coordinate system
+         * by: 
+         *       Bi_prime=[ dN/dx ]
+         *                [ dN/dy ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+        /*Build B': */
+        for(int i=0;i<numnodes;i++){
+                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
+                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
+
+}/*}}}*/
 void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
         element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
@@ -163 +305 @@
         xDelete<int>(doflist);
 }/*}}}*/
+
+/*Intermediaries*/
+void HydrologyShreveAnalysis::CreateHydrologyWaterVelocityInput(Element* element){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble dsdx,dsdy,dbdx,dbdy,w;
+
+        /*Retrieve all inputs and parameters*/
+        IssmDouble  rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+        IssmDouble  g         = element->GetMaterialParameter(ConstantsGEnum);
+        IssmDouble  CR        = element->GetMaterialParameter(HydrologyshreveCREnum);
+        IssmDouble  n_man     = element->GetMaterialParameter(HydrologyshreveNEnum);
+        IssmDouble  mu_water  = element->GetMaterialParameter(MaterialsMuWaterEnum);
+        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
+        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
+        Input* bedslopex_input     = element->GetInput(BedSlopeXEnum);     _assert_(bedslopex_input);
+        Input* bedslopey_input     = element->GetInput(BedSlopeYEnum);     _assert_(bedslopey_input);
+        Input* watercolumn_input   = element->GetInput(WatercolumnEnum);   _assert_(watercolumn_input);
+
+        /* compute VelocityFactor */
+        IssmDouble VelocityFactor = n_man*CR*CR*rho_water*g/mu_water;
+
+        /*Fetch number of vertices and allocate output*/
+        int numvertices = element->GetNumberOfVertices();
+        IssmDouble* vx  = xNew<IssmDouble>(numvertices);
+        IssmDouble* vy  = xNew<IssmDouble>(numvertices);
+
+        Gauss* gauss=element->NewGauss();
+        for(int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+                surfaceslopex_input->GetInputValue(&dsdx,gauss);
+                surfaceslopey_input->GetInputValue(&dsdy,gauss);
+                bedslopex_input->GetInputValue(&dbdx,gauss);
+                bedslopey_input->GetInputValue(&dbdy,gauss);
+                watercolumn_input->GetInputValue(&w,gauss);
+
+                /* Water velocity x and y components */
+        //      vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+        //      vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+                vx[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+                vy[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+        }
+
+        /*clean-up*/
+        delete gauss;
+
+        /*Add to inputs*/
+        element->AddInput(HydrologyWaterVxEnum,vx,P1Enum);
+        element->AddInput(HydrologyWaterVyEnum,vy,P1Enum);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(vy);
+}/*}}}*/

issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h ¶

@@ -25 +25 @@
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
+
+                /*Intermediaries*/
+                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+                void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+                void CreateHydrologyWaterVelocityInput(Element* element);
 };
 #endif

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

@@ -252 +252 @@
         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
-        IssmDouble     D[2][2];
+        IssmDouble     D[2][2]={0.};
 
         /*Retrieve all inputs and parameters*/
@@ -298 +298 @@
 
                 D[0][0]=D_scalar*dvxdx;
-                D[0][1]=0.;
-                D[1][0]=0.;
                 D[1][1]=D_scalar*dvydy;
                 TripleMultiply(B,2,numnodes,1,
@@ -326 +324 @@
                         vyaverage_input->GetInputAverage(&vy);
                         D[0][0]=h/2.0*fabs(vx);
-                        D[0][1]=0.;
-                        D[1][0]=0.;
                         D[1][1]=h/2.0*fabs(vy);
                 }

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp ¶

@@ -241 +241 @@
 
         switch(enum_in){
-                case MaterialsRhoIceEnum:                   return this->rho_ice;
-                case MaterialsRhoWaterEnum:                 return this->rho_water;
-                case MaterialsRhoFreshwaterEnum:            return this->rho_freshwater;
-                case MaterialsMuWaterEnum:                  return this->mu_water;
-                case MaterialsHeatcapacityEnum:             return this->heatcapacity;
-                case MaterialsThermalconductivityEnum:      return this->thermalconductivity;
-                case MaterialsTemperateiceconductivityEnum: return this->temperateiceconductivity;
-                case MaterialsLatentheatEnum:               return this->latentheat;
-                case MaterialsBetaEnum:                     return this->beta;
-                case MaterialsMeltingpointEnum:             return this->meltingpoint;
-                case ConstantsReferencetemperatureEnum:     return this->referencetemperature;
-                case MaterialsMixedLayerCapacityEnum:       return this->mixed_layer_capacity;
-                case MaterialsThermalExchangeVelocityEnum:  return this->thermal_exchange_velocity;
-                case HydrologydcSedimentPorosityEnum:       return this->sediment_porosity;
-                case HydrologydcSedimentThicknessEnum:      return this->sediment_thickness;
-                case HydrologydcSedimentCompressibilityEnum:return this->sediment_compressibility;
-                case HydrologydcEplPorosityEnum:            return this->epl_porosity;
-                case HydrologydcEplCompressibilityEnum:     return this->epl_compressibility;
-                case HydrologydcWaterCompressibilityEnum:   return this->water_compressibility;
-                case HydrologydcSedimentTransmitivityEnum:  return this->sediment_transmitivity;
-                case ConstantsGEnum:                        return this->g;
+                case MaterialsRhoIceEnum:                    return this->rho_ice;
+                case MaterialsRhoWaterEnum:                  return this->rho_water;
+                case MaterialsRhoFreshwaterEnum:             return this->rho_freshwater;
+                case MaterialsMuWaterEnum:                   return this->mu_water;
+                case MaterialsHeatcapacityEnum:              return this->heatcapacity;
+                case MaterialsThermalconductivityEnum:       return this->thermalconductivity;
+                case MaterialsTemperateiceconductivityEnum:  return this->temperateiceconductivity;
+                case MaterialsLatentheatEnum:                return this->latentheat;
+                case MaterialsBetaEnum:                      return this->beta;
+                case MaterialsMeltingpointEnum:              return this->meltingpoint;
+                case ConstantsReferencetemperatureEnum:      return this->referencetemperature;
+                case MaterialsMixedLayerCapacityEnum:        return this->mixed_layer_capacity;
+                case MaterialsThermalExchangeVelocityEnum:   return this->thermal_exchange_velocity;
+                case HydrologydcSedimentPorosityEnum:        return this->sediment_porosity;
+                case HydrologydcSedimentThicknessEnum:       return this->sediment_thickness;
+                case HydrologydcSedimentCompressibilityEnum: return this->sediment_compressibility;
+                case HydrologydcEplPorosityEnum:             return this->epl_porosity;
+                case HydrologydcEplCompressibilityEnum:      return this->epl_compressibility;
+                case HydrologydcWaterCompressibilityEnum:    return this->water_compressibility;
+                case HydrologydcSedimentTransmitivityEnum:   return this->sediment_transmitivity;
+                case HydrologyshreveCREnum:                  return this->hydro_CR;
+                case HydrologyshreveKnEnum:                  return this->hydro_kn;
+                case HydrologyshreveNEnum:                   return this->hydro_n;
+                case HydrologyshrevePEnum:                   return this->hydro_p;
+                case HydrologyshreveQEnum:                   return this->hydro_q;
+                case ConstantsGEnum:                         return this->g;
                 default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
         }