Changeset 20683

Timestamp:

06/01/16 17:05:50 (9 years ago)

Author:

felicity

Message:

CHG: added rigid body correction to vorticity in Budd et al. (2013) rheology

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Element.cpp (modified) (3 diffs)
• classes/Materials/Matearl.cpp (modified) (1 diff)
• shared/Elements/EarlComponents.cpp (modified) (3 diffs)
• shared/Elements/elements.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp ¶

@@ -224 +224 @@
 
         /*Intermediaries*/
-        IssmDouble vx,vy,vz;
-        IssmDouble dvx[3],dvy[3],dvz[3];
+        IssmDouble vx,vy,vz,vmag;
+        IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3];
         IssmDouble epso,epsprime; 
         int         dim;
@@ -241 +241 @@
         int numvertices = this->GetNumberOfVertices();
         IssmDouble* lambdas = xNew<IssmDouble>(numvertices);
+        IssmDouble* vorticityx = xNew<IssmDouble>(numvertices);
+        IssmDouble* vorticityy = xNew<IssmDouble>(numvertices);
+        IssmDouble* vorticityz = xNew<IssmDouble>(numvertices);
+        IssmDouble* omega_corrx = xNew<IssmDouble>(numvertices);
+        IssmDouble* omega_corry = xNew<IssmDouble>(numvertices);
+        IssmDouble* omega_corrz = xNew<IssmDouble>(numvertices);
 
         /* Start looping on the number of vertices: */
@@ -264 +270 @@
                         dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
                 }
-                EarlStrainrateQuantities(&epso,&epsprime,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+                /*Calculate velocity magnitude and its derivative*/
+                vmag = sqrt(vx*vx+vy*vy+vz*vz);
+                if(vmag<1e-12){
+                        vmag=1e-12;
+                        dvmag[0]=0;
+                        dvmag[1]=0;
+                        dvmag[2]=0;
+                }
+                else{
+                        dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
+                        dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
+                        dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
+                }
+                EarlStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
                 lambdas[iv]=EarlLambdaS(epso,epsprime);
+                vorticityx[iv]=dvz[1]-dvy[2];
+                vorticityy[iv]=dvx[2]-dvz[0];
+                vorticityz[iv]=dvy[0]-dvx[1];
+                omega_corrx[iv] = -vorticityx[iv] + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag;
+                omega_corry[iv] = -vorticityy[iv] + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag;
+                omega_corrz[iv] = -vorticityz[iv] + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag;
         }
 
         /*Add Stress tensor components into inputs*/
         this->AddInput(LambdaSEnum,lambdas,P1Enum);
-
+        this->AddInput(Outputdefinition21Enum,vorticityx,P1Enum);
+        this->AddInput(Outputdefinition22Enum,vorticityy,P1Enum);
+        this->AddInput(Outputdefinition23Enum,vorticityz,P1Enum);
+        this->AddInput(Outputdefinition31Enum,omega_corrx,P1Enum);
+        this->AddInput(Outputdefinition32Enum,omega_corry,P1Enum);
+        this->AddInput(Outputdefinition33Enum,omega_corrz,P1Enum);
+        
         /*Clean up and return*/
         delete gauss;

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

@@ -409 +409 @@
         IssmDouble E,lambdas;
         IssmDouble epso,epsprime_norm;
-        
-        EarlStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+        IssmDouble vmag,dvmag[3];
+
+        /*Calculate velocity magnitude and its derivative*/
+        vmag = sqrt(vx*vx+vy*vy+vz*vz);
+        if(vmag<1e-12){
+                dvmag[0]=0;
+                dvmag[1]=0;
+                dvmag[2]=0;
+        }
+        else{
+                dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
+                dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
+                dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
+        }
+
+        EarlStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
         lambdas=EarlLambdaS(epso,epsprime_norm);
 

issm/trunk-jpl/src/c/shared/Elements/EarlComponents.cpp ¶

@@ -2 +2 @@
 #include "../Numerics/types.h"
 #include "../Exceptions/exceptions.h"
-
-void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
+#include "./elements.h"
+void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
 
         /*Intermediaries*/
-        IssmDouble vorticity[3],vorticity_norm;
+        IssmDouble omega[3],omega_norm;
         IssmDouble nrsp[3],nrsp_norm;
         IssmDouble eps[3][3],epso;
         IssmDouble epsprime[3],epsprime_norm;
 
-        /*Create vorticity vector*/
-        _assert_(dvx && dvy && dvz);
-        vorticity[0] =  dvz[1] - dvy[2];
-        vorticity[1] =  dvx[2] - dvz[0];
-        vorticity[2] =  dvy[0] - dvx[1];
-
-        /*Normalize*/
-        vorticity_norm = sqrt(vorticity[0]*vorticity[0] + vorticity[1]*vorticity[1] + vorticity[2]*vorticity[2]);
-        if(vorticity_norm==0){
-                vorticity[0] = 0.;
-                vorticity[1] = 0.;
-                vorticity[2] = 1.;
-        }
-        else{
-                vorticity[0] =vorticity[0]/vorticity_norm;
-                vorticity[1] =vorticity[1]/vorticity_norm;
-                vorticity[2] =vorticity[2]/vorticity_norm;
-        }
+        /*Get omega, correction for rigid body rotation*/
+        EarlOmega(&omega[0],vx,vy,vz,vmag,dvx,dvy,dvz,dvmag);
 
         /*Non-rotating shear plane*/
-        nrsp[0] =  vy*vorticity[2] - vz*vorticity[1];
-        nrsp[1] =  vz*vorticity[0] - vx*vorticity[2];
-        nrsp[2] =  vx*vorticity[1] - vy*vorticity[0];
+        nrsp[0] =  vy*omega[2] - vz*omega[1];
+        nrsp[1] =  vz*omega[0] - vx*omega[2];
+        nrsp[2] =  vx*omega[1] - vy*omega[0];
 
         /*Normalize*/
@@ -81 +65 @@
                 for(int j=0;j<3;j++){
                         for(int k=0;k<3;k++){
-                                epsprime[j] += -nrsp[i]*eps[i][k]*vorticity[k]*vorticity[j];
+                                epsprime[j] += -nrsp[i]*eps[i][k]*omega[k]*omega[j];
                         }
                 }
@@ -88 +72 @@
         /*Get norm of epsprime*/
         epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
-
+        
         /*Assign output pointers*/
         *pepso=epso;
         *pepsprime_norm=epsprime_norm;
+}/*}}}*/
+void EarlOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble omega_norm;
+
+        /*Create vorticity vector*/
+        _assert_(dvx && dvy && dvz && dvmag);
+        if(vmag<1e-12)vmag=1e-12;
+
+        /*Create vorticity vector, corrected for rigid body rotation
+         * \overline{\omega} =\omega - \Omega
+         *                   =-\nabla\times{\bf v} - 2*\nabla U\times{\bf v}/U;
+         * check the magnitude of the second term -- if it is small, then the two
+         * vorticities (omega and first term in omega) are approx. equal
+         *
+         * */
+        omega[0] = -(dvz[1] - dvy[2]) + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag;
+        omega[1] = -(dvx[2] - dvz[0]) + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag;
+        omega[2] = -(dvy[0] - dvx[1]) + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag;
+
+        /*Take out vorticity component aligned with the velocity*/
+        IssmDouble wdotv = vx/vmag*omega[0] + vy/vmag*omega[1] + vz/vmag*omega[2];
+        omega[0] = omega[0] - wdotv*vx/vmag;
+        omega[1] = omega[1] - wdotv*vy/vmag;
+        omega[2] = omega[2] - wdotv*vz/vmag;
+
+        /*Normalize*/
+        omega_norm = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
+        if(omega_norm==0){
+                omega[0] = 0.;
+                omega[1] = 0.;
+                omega[2] = 1.;
+        }
+        else{
+                omega[0] =omega[0]/omega_norm;
+                omega[1] =omega[1]/omega_norm;
+                omega[2] =omega[2]/omega_norm;
+        }
+
 }/*}}}*/
 IssmDouble EarlLambdaS(IssmDouble epso, IssmDouble epsprime_norm){/*{{{*/

issm/trunk-jpl/src/c/shared/Elements/elements.h ¶

@@ -16 +16 @@
 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
 IssmDouble EarlLambdaS(IssmDouble epso, IssmDouble epsprime_norm);
-void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
+void EarlOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag);
+void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
 IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
                                  IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm,