Changeset 18496

Timestamp:

09/10/14 14:13:42 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added some forcings

File:

• issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp (modified) (8 diffs)

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

@@ -39 +39 @@
         iomodel->FetchDataToInput(elements,VxEnum);
         iomodel->FetchDataToInput(elements,VyEnum);
+        iomodel->FetchDataToInput(elements,SeaiceCoriolisFactorEnum);
 }/*}}}*/
 void SeaiceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -73 +74 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,D_scalar,dt;
+        IssmDouble  Jdet,D_scalar,dt,thickness,concentration;
         IssmDouble  C[3][3];
         IssmDouble* xyz_list = NULL;
@@ -89 +90 @@
         element->GetVerticesCoordinates(&xyz_list);
         element->FindParam(&dt,TimesteppingTimeStepEnum);
+        IssmDouble rho_i               = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        Input*     thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -95 +98 @@
 
                 gauss->GaussPoint(ig);
+
+                /*Get Jacobian Determinant and thickness*/
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                thickness_input->GetInputValue(&thickness,gauss);
+                concentration_input->GetInputValue(&concentration,gauss);
 
                 /*Create first part of the stiffness matrix*/
                 this->GetM(M,element,xyz_list,gauss);
-                D_scalar=gauss->weight*Jdet;
+                D_scalar=rho_i*thickness*gauss->weight*Jdet;
                 TripleMultiply(M,1,numdof,1,
                                         &D_scalar,1,1,0,
@@ -117 +124 @@
                                         B,3,numdof,0,
                                         &Ke->values[0],1);
+
+                /*Create Coriolis part*/
         }
 
@@ -133 +142 @@
 
         /*Intermediaries */
-        IssmDouble  concentration,thickness,dt,Jdet,vx,vy;
-        IssmDouble  rho_air,air_coef,air_lin_drag_coef,air_quad_drag_coef;
-        IssmDouble  windvx,windvy,wind_vnorm;
+        IssmDouble  air_coef,ocean_coef,constant_part;
+        IssmDouble  rho_ice,rho_air,rho_ocean,gravity;
+        IssmDouble  vx,vy,vxstar,vystar,windvx,windvy,oceanvx,oceanvy,vnorm;
+        IssmDouble  air_lin_drag_coef,air_quad_drag_coef;
+        IssmDouble  ocean_lin_drag_coef,ocean_quad_drag_coef;
+        IssmDouble  concentration,thickness,coriolis_factor,dt,Jdet;
+        IssmDouble  ocean_turning_angle,dssh[2];
         IssmDouble* xyz_list = NULL;
 
@@ -153 +166 @@
         element->FindParam(&air_lin_drag_coef,SurfaceforcingsAirLinDragCoefEnum);
         element->FindParam(&air_quad_drag_coef,SurfaceforcingsAirQuadDragCoefEnum);
+        element->FindParam(&rho_ocean,BasalforcingsRhoOceanEnum);
+        element->FindParam(&ocean_coef,BasalforcingsOceanCoefEnum);
+        element->FindParam(&ocean_lin_drag_coef,BasalforcingsOceanLinDragCoefEnum);
+        element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
+        element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
+        element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
+        rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        gravity = element->GetMaterialParameter(ConstantsGEnum);
+        Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
+        Input* coriolis_fact_input = element->GetInput(SeaiceCoriolisFactorEnum);  _assert_(coriolis_fact_input);
         Input* vx_input            = element->GetInput(VxEnum);                    _assert_(vx_input);
         Input* vy_input            = element->GetInput(VyEnum);                    _assert_(vy_input);
+        Input* vxstar_input        = element->GetInput(VxStarEnum);                _assert_(vxstar_input);
+        Input* vystar_input        = element->GetInput(VyStarEnum);                _assert_(vystar_input);
         Input* windvx_input        = element->GetInput(SurfaceforcingsWindVxEnum); _assert_(windvx_input);
         Input* windvy_input        = element->GetInput(SurfaceforcingsWindVyEnum); _assert_(windvy_input);
-        Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
+        Input* oceanvx_input       = element->GetInput(BasalforcingsOceanVxEnum);  _assert_(oceanvx_input);
+        Input* oceanvy_input       = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
+        Input* oceanssh_input      = element->GetInput(BasalforcingsOceanSshEnum); _assert_(oceanssh_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -168 +195 @@
                 element->NodalFunctions(basis, gauss);
 
+                /*Get all inputs on gauss point*/
                 thickness_input->GetInputValue(&thickness,gauss);
                 concentration_input->GetInputValue(&concentration,gauss);
+                coriolis_fact_input->GetInputValue(&coriolis_factor,gauss);
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
+                vxstar_input->GetInputValue(&vxstar,gauss);
+                vystar_input->GetInputValue(&vystar,gauss);
                 windvx_input->GetInputValue(&windvx,gauss);
                 windvy_input->GetInputValue(&windvy,gauss);
+                oceanvx_input->GetInputValue(&oceanvx,gauss);
+                oceanvy_input->GetInputValue(&oceanvy,gauss);
+                oceanssh_input->GetInputDerivativeValue(&dssh[2],xyz_list,gauss);
+
+                /*Previous speed (finite differencing)*/
+                for(int i=0;i<numnodes;i++){
+                        pe->values[i*2+0]+=Jdet*gauss->weight*rho_ice*thickness*vx*basis[i];
+                        pe->values[i*2+1]+=Jdet*gauss->weight*rho_ice*thickness*vy*basis[i];
+                }
 
                 /*Atmespheric forcing*/
-                wind_vnorm = sqrt(windvx*windvx + windvy*windvy);
-
+                vnorm    = sqrt(windvx*windvx + windvy*windvy);
+                constant_part = concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*vnorm);
                 for(int i=0;i<numnodes;i++){
-                        pe->values[i*2+0] += dt*concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*wind_vnorm)*windvx*Jdet*gauss->weight*basis[i];
-                        pe->values[i*2+1] += dt*concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*wind_vnorm)*windvy*Jdet*gauss->weight*basis[i];
+                        pe->values[i*2+0] += dt*constant_part*windvx*Jdet*gauss->weight*basis[i];
+                        pe->values[i*2+1] += dt*constant_part*windvy*Jdet*gauss->weight*basis[i];
+                }
+
+                /*Ocean forcing*/
+                vnorm = sqrt(pow(oceanvx-vx,2) + pow(oceanvy-vy,2));
+                constant_part = concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*vnorm);
+                for(int i=0;i<numnodes;i++){
+                        pe->values[i*2+0] += dt*constant_part*(oceanvy-vy)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+                        pe->values[i*2+1] += dt*constant_part*(vx-oceanvx)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+
+                        pe->values[i*2+0] += dt*constant_part*oceanvx*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+                        pe->values[i*2+1] += dt*constant_part*oceanvy*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+                }
+
+                /*Coriolis forces (use ustar)*/
+                for(int i=0;i<numnodes;i++){
+                        pe->values[i*2+0] += dt*(-rho_ice*thickness*coriolis_factor*vystar)*Jdet*gauss->weight*basis[i];
+                        pe->values[i*2+1] += dt*(+rho_ice*thickness*coriolis_factor*vxstar)*Jdet*gauss->weight*basis[i];
+
+                        pe->values[i*2+0] += dt*(-rho_ice*thickness*gravity*dssh[0])*Jdet*gauss->weight*basis[i];
+                        pe->values[i*2+1] += dt*(-rho_ice*thickness*gravity*dssh[1])*Jdet*gauss->weight*basis[i];
                 }
         }