Changeset 22264

Timestamp:

11/16/17 14:05:37 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: trying out a new method to interpolate ice thickness based on diffusion

File:

• issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (modified) (7 diffs)

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

@@ -71 +71 @@
 
         /*Intermediaries */
-        int         n=3;
-        IssmDouble  Jdet,omega,ds[2],slope;
+        IssmDouble  yts = 365*24*3600.;
+        IssmDouble  Jdet,vx,vy,vel;
         IssmDouble* xyz_list = NULL;
 
@@ -84 +84 @@
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-        Input* omega_input         = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
-        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); 
-        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); 
-        IssmDouble rhog    = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+        Input* vx_input = element->GetInput(SurfaceSlopeXEnum); _assert_(vx_input); 
+        Input* vy_input = element->GetInput(SurfaceSlopeYEnum); _assert_(vy_input);
+
+        /*make sure are diffusivisty is large enough*/
+        vel = sqrt(vx*vx+vy*vy);
+        if(sqrt(vx*vx+vy*vy)==0.){
+                vx = 0.1/yts;
+                vy = 0.1/yts;
+        }
+        else if(vel<0.1/yts){
+                vx = vx/vel*0.1;
+                vy = vy/vel*0.1;
+
+        }
 
         /* Start  looping on the number of gaussian points: */
@@ -95 +105 @@
                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                omega_input->GetInputValue(&omega,gauss);
-                surfaceslopex_input->GetInputValue(&ds[0],gauss);
-                surfaceslopey_input->GetInputValue(&ds[1],gauss);
-
-                slope = sqrt(ds[0]*ds[0] + ds[1]*ds[1]);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
 
                 for(int i=0;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
-                                Ke->values[i*numnodes+j] += pow(rhog,n)*exp(omega)*pow(slope,n-1)*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(vx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
                         }
                 }
@@ -117 +124 @@
 
         /*Intermediaries */
-        IssmDouble  dhdt,mb,ms,Jdet;
+        IssmDouble  dhdt[2],mb[2],ms[2],Jdet;
         IssmDouble* xyz_list = NULL;
 
@@ -141 +148 @@
                 element->NodalFunctions(basis,gauss);
 
-                ms_input->GetInputValue(&ms,gauss);
-                mb_input->GetInputValue(&mb,gauss);
-                dhdt_input->GetInputValue(&dhdt,gauss);
+                ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
+
+                ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
+                mb_input->GetInputDerivativeValue(&mb[0],xyz_list,gauss);
+                dhdt_input->GetInputDerivativeValue(&dhdt[0],xyz_list,gauss);
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
-                                        (ms-mb-dhdt)*basis[i]
+                                        (ms[0]+ms[1]-mb[0]-mb[1]-dhdt[0]-dhdt[1])*basis[i]
                                         );
         }
@@ -157 +166 @@
 }/*}}}*/
 void           Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-                element->GetSolutionFromInputsOneDof(solution,SurfaceEnum);
+                element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
 }/*}}}*/
 void           Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
@@ -164 +173 @@
 void           Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
-        /*Intermediaries*/
-        IssmDouble  ds[2],s,b,D;
-        IssmDouble* xyz_list = NULL;
+                        element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
 
-        //element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
-        element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
-
-        /*Fetch number of vertices and allocate velocity vectors*/
-        int numvertices = element->GetNumberOfVertices();
-        IssmDouble* vel_list = xNew<IssmDouble>(numvertices);
-        IssmDouble* vx_list  = xNew<IssmDouble>(numvertices);
-        IssmDouble* vy_list  = xNew<IssmDouble>(numvertices);
-
-        /*Retrieve all inputs and parameters*/
-        element->GetVerticesCoordinates(&xyz_list);
-        Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);
-        Input* H_input   = element->GetInput(ThicknessEnum);                            _assert_(H_input);
-        Input* s_input   = element->GetInput(SurfaceEnum);                              _assert_(s_input);
-        Input* b_input   = element->GetInput(BaseEnum);                                 _assert_(b_input);
-
-        /*Calculate velocities*/
-        Gauss* gauss=element->NewGauss();
-        for(int iv=0;iv<numvertices;iv++){
-                gauss->GaussVertex(iv);
-
-                if(D_input){
-                        D_input->GetInputValue(&D,gauss);
-                }
-                else{
-                        D = 0.;
-                }
-                b_input->GetInputValue(&b,gauss);
-                s_input->GetInputValue(&s,gauss);
-                s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
-
-                vx_list[iv] = -1./(s-b)*D*ds[0];
-                vy_list[iv] = -1./(s-b)*D*ds[1];
-                vel_list[iv] = sqrt(pow(vx_list[iv],2) + pow(vy_list[iv],2));
-        }
-
-        /*Add vx and vy as inputs to the tria element: */
-        element->AddInput(VxEnum,vx_list,P1Enum);
-        element->AddInput(VyEnum,vy_list,P1Enum);
-        element->AddInput(VelEnum,vel_list,P1Enum);
-
-        /*Free ressources:*/
-        delete gauss;
-        xDelete<IssmDouble>(vy_list);
-        xDelete<IssmDouble>(vx_list);
-        xDelete<IssmDouble>(vel_list);
-        xDelete<IssmDouble>(xyz_list);
 }/*}}}*/
 void           Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/