Changeset 11361

Timestamp:

02/08/12 13:29:22 (13 years ago)

Author:

seroussi

Message:

fixing stabilization for enthalpy solution

Location:

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

Files:

• Penta.cpp (modified) (14 diffs)
• Penta.h (modified) (1 diff)

issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

@@ -1132 +1132 @@
 /*}}}*/
 /*FUNCTION Penta::GetStabilizationParameter {{{1*/
-double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
+double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double kappa){
         /*Compute stabilization parameter*/
+        /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
+        /*kappa=enthalpydiffusionparameter for enthalpy model*/
 
         double normu;
@@ -1139 +1141 @@
 
         normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
-        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
-                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
+        if(normu*diameter/(3*2*kappa)<1){ 
+                tau_parameter=pow(diameter,2)/(3*2*2*kappa);
         }
         else tau_parameter=diameter/(2*normu);
@@ -3372 +3374 @@
                 else if(stabilization==2){
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
-                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
+                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
 
                         for(i=0;i<numdof;i++){
@@ -3486 +3488 @@
         double     Jdet,u,v,w,um,vm,wm,vel;
         double     h,hx,hy,hz,vx,vy,vz;
-        double     gravity,rho_ice,rho_water;
+        double     gravity,rho_ice,rho_water,kappa;
         double     heatcapacity,thermalconductivity,dt;
         double     tau_parameter,diameter;
@@ -3512 +3514 @@
         heatcapacity=matpar->GetHeatCapacity();
         thermalconductivity=matpar->GetThermalConductivity();
+        kappa=thermalconductivity/(rho_ice*heatcapacity);
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
@@ -3603 +3606 @@
                 else if(stabilization==2){
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
-                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
+                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
 
                         for(i=0;i<numdof;i++){
@@ -3708 +3711 @@
         double Jdet,phi,dt;
         double rho_ice,heatcapacity;
-        double thermalconductivity;
-        double viscosity,enthalpy;
+        double thermalconductivity,kappa;
+        double viscosity,enthalpy,pressure;
         double tau_parameter,diameter;
         double u,v,w;
@@ -3733 +3736 @@
         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
+        Input* pressure_input=inputs->GetInput(PressureEnum);      _assert_(pressure_input);
         Input* enthalpy_input=NULL;
         if (dt) enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(inputs);
@@ -3746 +3750 @@
                 GetNodalFunctionsP1(&L[0], gauss);
 
+                enthalpy_input->GetInputValue(&enthalpy, gauss);
                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
                 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
@@ -3757 +3762 @@
                 /* Build transient now */
                 if(dt){
-                        enthalpy_input->GetInputValue(&enthalpy, gauss);
                         scalar_transient=enthalpy*Jdet*gauss->weight;
                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
@@ -3768 +3772 @@
                         vy_input->GetInputValue(&v, gauss);
                         vz_input->GetInputValue(&w, gauss);
-
-                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+                        pressure_input->GetInputValue(&pressure, gauss);
+                        kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
+
+                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
 
                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
@@ -3961 +3967 @@
         double Jdet,phi,dt;
         double rho_ice,heatcapacity;
-        double thermalconductivity;
+        double thermalconductivity,kappa;
         double viscosity,temperature;
         double tau_parameter,diameter;
@@ -3981 +3987 @@
         heatcapacity=matpar->GetHeatCapacity();
         thermalconductivity=matpar->GetThermalConductivity();
+        kappa=heatcapacity/(rho_ice*thermalconductivity);
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
@@ -4022 +4029 @@
                         vz_input->GetInputValue(&w, gauss);
 
-                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
 
                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);

issm/trunk-jpl/src/c/objects/Elements/Penta.h

@@ -183 +183 @@
                 void      GetPhi(double* phi, double*  epsilon, double viscosity);
                 void      GetSolutionFromInputsEnthalpy(Vec solutiong);
-                double  GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity);
+                double  GetStabilizationParameter(double u, double v, double w, double diameter, double kappa);
                 void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
                 void    GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);