Changeset 16895

Timestamp:

11/22/13 13:58:57 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixed Enthalpy

Location:

issm/trunk-jpl/src/c/analyses

Files:

• EnthalpyAnalysis.cpp (modified) (10 diffs)
• EnthalpyAnalysis.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -236 +236 @@
 
         /*Enthalpy diffusion parameter*/
-        IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element); _assert_(kappa>0.);
+        IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyEnum); _assert_(kappa>0.);
 
         /* Start  looping on the number of gaussian points: */
@@ -304 +304 @@
                                 for(int j=0;j<numnodes;j++){
                                         Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
-                                          ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]);
+                                          ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
                                 }
                         }
@@ -310 +310 @@
                                 for(int i=0;i<numnodes;i++){
                                         for(int j=0;j<numnodes;j++){
-                                                Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]);
+                                                Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
                                         }
                                 }
@@ -410 +410 @@
         IssmDouble*    basis          = xNew<IssmDouble>(numnodes);
         IssmDouble*    dbasis         = xNew<IssmDouble>(3*numnodes);
-        IssmDouble*    pressure       = xNew<IssmDouble>(numvertices);
-        IssmDouble*    enthalpypicard = xNew<IssmDouble>(numvertices);
 
         /*Retrieve all inputs and parameters*/
@@ -425 +423 @@
         if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
         if(stabilization==2){
-                element->GetInputListOnVertices(enthalpypicard,EnthalpyPicardEnum);
-                element->GetInputListOnVertices(pressure,PressureEnum);
                 diameter=element->MinEdgeLength(xyz_list);
+                kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.);
         }
 
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=element->NewGauss(3);
+        Gauss* gauss=element->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
@@ -440 +437 @@
 
                 scalar_def=phi/rho_ice*Jdet*gauss->weight;
-                if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
+                if(dt!=0.) scalar_def=scalar_def*dt;
 
                 /*TODO: add -beta*laplace T_m(p)*/
@@ -458 +455 @@
                         vy_input->GetInputValue(&v,gauss);
                         vz_input->GetInputValue(&w,gauss);
-                        kappa          = element->EnthalpyDiffusionParameterVolume(numvertices,enthalpypicard,pressure) / rho_ice;
-                        tau_parameter  = element->StabilizationParameter(u,v,w,diameter,kappa);
-
-                        for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
-                        if(reCast<bool,IssmDouble>(dt)){
-                                for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+                        tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
+
+                        for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
+
+                        if(dt!=0.){
+                                for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
                         }
                 }
@@ -471 +468 @@
         xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(dbasis);
-        xDelete<IssmDouble>(pressure);
-        xDelete<IssmDouble>(enthalpypicard);
         xDelete<IssmDouble>(xyz_list);
         delete gauss;
@@ -798 +793 @@
         }
 }/*}}}*/
-IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element){/*{{{*/
+IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
 
         int         iv;
@@ -813 +808 @@
         IssmDouble* dHpmp       = xNew<IssmDouble>(numvertices);
         element->GetInputListOnVertices(pressures,PressureEnum);
-        element->GetInputListOnVertices(enthalpies,EnthalpyEnum);
+        element->GetInputListOnVertices(enthalpies,enthalpy_enum);
         for(iv=0;iv<numvertices;iv++){
                 PIE[iv]   = PureIceEnthalpy(element,pressures[iv]);

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h

@@ -36 +36 @@
                 /*Intermediaries*/
                 IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
-                IssmDouble EnthalpyDiffusionParameterVolume(Element* element);
+                IssmDouble EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum);
                 IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
                 IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);