Changeset 18501

Timestamp:

09/10/14 21:06:58 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added ocean forcing to stiffness matrix

File:

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

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

@@ -75 +75 @@
         /*Intermediaries */
         IssmDouble  Jdet,D_scalar,dt,thickness,concentration;
+        IssmDouble  ocean_coef,rho_ocean,ocean_lin_drag_coef,ocean_quad_drag_coef;
+        IssmDouble  vx,vy,oceanvx,oceanvy,vnorm,ocean_turning_angle;
         IssmDouble  C[3][3];
         IssmDouble* xyz_list = NULL;
@@ -90 +92 @@
         element->GetVerticesCoordinates(&xyz_list);
         element->FindParam(&dt,TimesteppingTimeStepEnum);
-        IssmDouble rho_i               = element->GetMaterialParameter(MaterialsRhoIceEnum);
-        Input*     thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
+        element->FindParam(&ocean_coef,BasalforcingsOceanCoefEnum);
+        element->FindParam(&ocean_lin_drag_coef,BasalforcingsOceanLinDragCoefEnum);
+        element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
+        element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
+        IssmDouble rho_i           = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        Input*     thickness_input = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
+        Input*     vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
+        Input*     vy_input        = element->GetInput(VyEnum);                    _assert_(vy_input);
+        Input*     oceanvx_input   = element->GetInput(BasalforcingsOceanVxEnum);  _assert_(oceanvx_input);
+        Input*     oceanvy_input   = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -125 +135 @@
                                         &Ke->values[0],1);
 
-                /*Create Coriolis part*/
+                /*Create Ocean forcing part*/
+                oceanvx_input->GetInputValue(&oceanvx,gauss);
+                oceanvy_input->GetInputValue(&oceanvy,gauss);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                vnorm = sqrt(pow(oceanvx-vx,2) + pow(oceanvy-vy,2));
+                D_scalar = dt*concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm)*cos(ocean_turning_angle)*gauss->weight*Jdet;
+                TripleMultiply(M,1,numdof,1,
+                                        &D_scalar,1,1,0,
+                                        M,1,numdof,0,
+                                        &Ke->values[0],1);
         }
 
@@ -170 +190 @@
         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);
@@ -226 +245 @@
                 /*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);
+                constant_part = concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_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];