Changeset 17976 for issm/trunk-jpl/src/c/classes/FemModel.cpp

Timestamp:

05/09/14 15:37:18 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: pushing misfits away from elements

File:

• issm/trunk-jpl/src/c/classes/FemModel.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1134 +1134 @@
 void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/
 
-        IssmDouble J = 0;
+        /*output: */
+        IssmDouble J=0.;
         IssmDouble J_sum;
 
-        for(int i=0;i<this->elements->Size();i++){
-                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-                J+=element->BalancethicknessMisfit();
-        }
+        IssmDouble  weight,vx,vy,H,dvx[2],dvy[2],dH[2];
+        IssmDouble  temp,Jdet,dhdt,basal_melting,surface_mass_balance;
+        IssmDouble* xyz_list = NULL;
+        IssmDouble  dp[3];
+
+        /*Compute Misfit: */
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+
+                /*If on water, return 0: */
+                if(!element->IsIceInElement()) continue;
+
+                /* Get node coordinates*/
+                element->GetVerticesCoordinates(&xyz_list);
+                Input* weights_input              = element->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
+                Input* thickness_input            = element->GetInput(ThicknessEnum); _assert_(thickness_input);
+                Input* vx_input                   = element->GetInput(VxEnum);                                 _assert_(vx_input);
+                Input* vy_input                   = element->GetInput(VyEnum);                                 _assert_(vy_input);
+                Input* surface_mass_balance_input = element->GetInput(SurfaceforcingsMassBalanceEnum);         _assert_(surface_mass_balance_input);
+                Input* basal_melting_input        = element->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
+                Input* dhdt_input                 = element->GetInput(BalancethicknessThickeningRateEnum);     _assert_(dhdt_input);
+
+                /* Start  looping on the number of gaussian points: */
+                Gauss* gauss=element->NewGauss(2);
+                for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                        gauss->GaussPoint(ig);
+
+                        /* Get Jacobian determinant: */
+                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+                        /*Get all parameters at gaussian point*/
+                        weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum);
+                        thickness_input->GetInputValue(&H, gauss);
+                        thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
+                        surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
+                        basal_melting_input->GetInputValue(&basal_melting,gauss);
+                        dhdt_input->GetInputValue(&dhdt,gauss);
+                        vx_input->GetInputValue(&vx,gauss);
+                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                        vy_input->GetInputValue(&vy,gauss);
+                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+
+                        /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/
+                        temp  = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-basal_melting-dhdt);
+                        J    +=weight*1/2*temp*temp*Jdet*gauss->weight;
+                }
+
+                /*clean up and Return: */
+                xDelete<IssmDouble>(xyz_list);
+                delete gauss;
+        }
+
+        /*Sum all J from all cpus of the cluster:*/
         ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
@@ -1151 +1202 @@
 void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/
 
-        /*Intermediary*/
-        int i;
-        Element* element=NULL;
-
         /*output: */
-        IssmDouble J=0;
+        IssmDouble J=0.;
         IssmDouble J_sum;
 
+        IssmDouble  thickness,thicknessobs,weight;
+        IssmDouble  Jdet;
+        IssmDouble* xyz_list = NULL;
+        IssmDouble  dp[3];
+
         /*Compute Misfit: */
-        for (i=0;i<elements->Size();i++){
-                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                J+=element->ThicknessAbsGradient();
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+
+                /*If on water, return 0: */
+                if(!element->IsIceInElement()) continue;
+
+                /* Get node coordinates*/
+                element->GetVerticesCoordinates(&xyz_list);
+
+                /*Retrieve all inputs we will be needing: */
+                Input* weights_input   =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+                Input* thickness_input =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
+
+                /* Start  looping on the number of gaussian points: */
+                Gauss* gauss=element->NewGauss(2);
+                for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                        gauss->GaussPoint(ig);
+
+                        /* Get Jacobian determinant: */
+                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+                        /*Get all parameters at gaussian point*/
+                        weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum);
+                        thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
+
+                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+                        J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
+                }
+
+                /*clean up and Return: */
+                xDelete<IssmDouble>(xyz_list);
+                delete gauss;
         }