Changeset 14660

Timestamp:

04/19/13 13:48:09 (12 years ago)

Author:

bdef

Message:

CHG: Adding parameters to the HydroDC matrix construction

Location:

issm/trunk-jpl/src/c/classes/objects

Files:

• Elements/Tria.cpp (modified) (8 diffs)
• Materials/Matpar.h (modified) (1 diff)

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

@@ -5758 +5758 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/
+/*FUNCTION Tria::CreateKMatrixHydrologyShreve(void){{{*/
 ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){
-
+        
         /*Constants*/
         const int  numdof=NDOF1*NUMVERTICES;
-
-        /*Intermediaries */
+        
+/*Intermediaries */
         IssmDouble diffusivity;
         IssmDouble Jdettria,DL_scalar,dt,h;
@@ -5779 +5779 @@
         IssmDouble DLprime[2][2]                   ={0.0};
         GaussTria *gauss=NULL;
-
-        /*Skip if water or ice shelf element*/
+        
+/*Skip if water or ice shelf element*/
         if(IsOnWater() | IsFloating()) return NULL;
-
-        /*Initialize Element matrix*/
+        
+/*Initialize Element matrix*/
         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
-        /*Create water velocity vx and vy from current inputs*/
+        
+/*Create water velocity vx and vy from current inputs*/
         CreateHydrologyWaterVelocityInput();
-
-        /*Retrieve all inputs and parameters*/
+        
+/*Retrieve all inputs and parameters*/
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
@@ -5796 +5796 @@
         Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
         h=sqrt(2*this->GetArea());
-
-        /* Start  looping on the number of gaussian points: */
+        
+/* Start  looping on the number of gaussian points: */
         gauss=new GaussTria(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
-
+                
                 gauss->GaussPoint(ig);
-
+                
                 GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
                 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
+                
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
                 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
                 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
-
+                
                 DL_scalar=gauss->weight*Jdettria;
-
+                
                 TripleMultiply( &L[0],1,numdof,1,
-                                        &DL_scalar,1,1,0,
-                                        &L[0],1,numdof,0,
-                                        &Ke->values[0],1);
-
+                                                                                &DL_scalar,1,1,0,
+                                                                                &L[0],1,numdof,0,
+                                                                                &Ke->values[0],1);
+                
                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
-
+                
                 dvxdx=dvx[0];
                 dvydy=dvy[1];
                 DL_scalar=dt*gauss->weight*Jdettria;
-
+                
                 DL[0][0]=DL_scalar*dvxdx;
                 DL[1][1]=DL_scalar*dvydy;
                 DLprime[0][0]=DL_scalar*vx;
                 DLprime[1][1]=DL_scalar*vy;
-
+                
                 TripleMultiply( &B[0][0],2,numdof,1,
-                                        &DL[0][0],2,2,0,
-                                        &B[0][0],2,numdof,0,
-                                        &Ke->values[0],1);
-
+                                                                                &DL[0][0],2,2,0,
+                                                                                &B[0][0],2,numdof,0,
+                                                                                &Ke->values[0],1);
+                
                 TripleMultiply( &B[0][0],2,numdof,1,
-                                        &DLprime[0][0],2,2,0,
-                                        &Bprime[0][0],2,numdof,0,
-                                        &Ke->values[0],1);
-
+                                                                                &DLprime[0][0],2,2,0,
+                                                                                &Bprime[0][0],2,numdof,0,
+                                                                                &Ke->values[0],1);
+                
                 /*Artificial diffusivity*/
                 vel=sqrt(vx*vx+vy*vy);
@@ -5850 +5850 @@
                 KDL[0][1]=DL_scalar*K[0][1];
                 KDL[1][1]=DL_scalar*K[1][1];
-
+                
                 TripleMultiply( &Bprime[0][0],2,numdof,1,
-                                        &KDL[0][0],2,2,0,
-                                        &Bprime[0][0],2,numdof,0,
-                                        &Ke->values[0],1);
-        }
-
-        /*Clean up and return*/
+                                                                                &KDL[0][0],2,2,0,
+                                                                                &Bprime[0][0],2,numdof,0,
+                                                                                &Ke->values[0],1);
+        }
+        
+/*Clean up and return*/
         delete gauss;
         return Ke;
@@ -5870 +5870 @@
         /* Intermediaries */
         IssmDouble  D_scalar,Jdet;
-        IssmDouble  sediment_porosity,dt;
+        IssmDouble  sediment_compressibility, sediment_porosity, sediment_thickness,
+                sediment_transmitivity, water_compressibility, rho_freshwater, gravity, dt;
+        IssmDouble  sediment_storing;
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  B[2][numdof];
-        IssmDouble L[NUMVERTICES];
+        IssmDouble  L[NUMVERTICES];
         IssmDouble  D[2][2];
-        GaussTria  *gauss = NULL;
+        GaussTria   *gauss = NULL;
 
         /*Initialize Element matrix*/
@@ -5882 +5884 @@
         /*Retrieve all inputs and parameters*/
         GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
-        sediment_porosity = matpar->GetSedimentPorosity();
+        sediment_compressibility                        = matpar->GetSedimentCompressibility();
+        sediment_porosity                                                       = matpar->GetSedimentPorosity();
+        sediment_thickness                                              = matpar->GetSedimentThickness();
+        sediment_transmitivity                          = matpar->GetSedimentTransmitivity();
+        water_compressibility                                   = matpar->GetWaterCompressibility();
+        rho_freshwater                                                          = matpar->GetRhoFreshwater();
+        gravity                                                                                         = matpar->GetG();
+        
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+        /*Compute Storing coefficient fro the above */
+
+        sediment_storing = rho_freshwater*gravity*sediment_porosity*sediment_thickness*
+                (water_compressibility+(sediment_compressibility/sediment_porosity));
 
         /* Start looping on the number of gaussian points: */
@@ -5900 +5914 @@
                 GetBHydro(&B[0][0],&xyz_list[0][0],gauss); 
                 TripleMultiply(&B[0][0],2,numdof,1,
-                                        &D[0][0],2,2,0,
-                                        &B[0][0],2,numdof,0,
-                                        &Ke->values[0],1);
+                                                                         &D[0][0],2,2,0,
+                                                                         &B[0][0],2,numdof,0,
+                                                                         &Ke->values[0],1);
 
                 /*Transient*/
@@ -5908 +5922 @@
                         GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
                         D_scalar=gauss->weight*Jdet;
-
+                        
                         TripleMultiply(&L[0],numdof,1,0,
-                                                &D_scalar,1,1,0,
-                                                &L[0],1,numdof,0,
-                                                &Ke->values[0],1);
+                                                                                 &D_scalar,1,1,0,
+                                                                                 &L[0],1,numdof,0,
+                                                                                 &Ke->values[0],1);
                 }
         }

issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h

@@ -111 +111 @@
                 IssmDouble GetSedimentTransmitivity();
                 IssmDouble GetWaterCompressibility();
-                IssmDouble GetWaterDensity();
                 IssmDouble TMeltingPoint(IssmDouble pressure);
                 IssmDouble PureIceEnthalpy(IssmDouble pressure);