Index: /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp (revision 18495) +++ /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp (revision 18496) @@ -39,4 +39,5 @@ iomodel->FetchDataToInput(elements,VxEnum); iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,SeaiceCoriolisFactorEnum); }/*}}}*/ void SeaiceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ @@ -73,5 +74,5 @@ /*Intermediaries */ - IssmDouble Jdet,D_scalar,dt; + IssmDouble Jdet,D_scalar,dt,thickness,concentration; IssmDouble C[3][3]; IssmDouble* xyz_list = NULL; @@ -89,4 +90,6 @@ 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,9 +98,13 @@ 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,4 +124,6 @@ B,3,numdof,0, &Ke->values[0],1); + + /*Create Coriolis part*/ } @@ -133,7 +142,11 @@ /*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,9 +166,23 @@ 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,17 +195,50 @@ 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;ivalues[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;ivalues[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;ivalues[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;ivalues[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]; } }