Changeset 18393

Timestamp:

08/14/14 16:15:31 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added dJ/dnux for the velocities

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointBalancethickness2Analysis.cpp (modified) (6 diffs)
• classes/FemModel.cpp (modified) (2 diffs)

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

@@ -234 +234 @@
         /*Intermediaries*/
         IssmDouble Jdet,weight;
-        IssmDouble vxobs,vyobs;
-        IssmDouble nux,nuy,dphi[2];
+        IssmDouble vxobs,vyobs,thicknessobs;
+        IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
         IssmDouble *xyz_list= NULL;
 
@@ -255 +255 @@
         Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
         Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
+        Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
 
         Gauss* gauss=element->NewGauss(2);
@@ -268 +269 @@
                 nuy_input->GetInputValue(&nuy,gauss);
                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+                thicknessobs_input->GetInputValue(&thicknessobs,gauss);
+
+                dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
+                nuvobs2   = nux*vxobs*nux*vxobs + nuy*vyobs*nuy*vyobs;
 
                 /*Build gradient vector (actually -dJ/da): */
                 for(int i=0;i<numvertices;i++){
-                        ge[i]+= - weight*Jdet*gauss->weight*(-vxobs)*dphi[1]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
+                        //ge[i]+= - weight*Jdet*gauss->weight*(-vxobs)*dphi[1]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i]; //THIS IS FOR Q//V
+                        ge[i]+= - weight*Jdet*gauss->weight*basis[i]*(
+                                                -2.*nux*vxobs*vxobs*(
+                                                        dphinorm2*dphinorm2 - thicknessobs*thicknessobs*dphinorm2*nuvobs2
+                                                        )/ ( (pow(nux,6)+pow(nuy,6))*pow(vxobs,6) + 3*(nux*vxobs*nuy*vyobs)*(nux*vxobs*nuy*vyobs)*nuvobs2)
+                                                );
                         _assert_(!xIsNan<IssmDouble>(ge[i]));
                 }
@@ -291 +301 @@
         /*Intermediaries*/
         IssmDouble Jdet,weight;
-        IssmDouble vxobs,vyobs;
-        IssmDouble nux,nuy,dphi[2];
+        IssmDouble vxobs,vyobs,thicknessobs;
+        IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
         IssmDouble *xyz_list= NULL;
 
@@ -306 +316 @@
         element->GetVerticesCoordinates(&xyz_list);
         element->GradientIndexing(&vertexpidlist[0],control_index);
-        Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
-        Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
-        Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
-        Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
-        Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
-        Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
+        Input* weights_input      = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+        Input* vxobs_input        = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
+        Input* vyobs_input        = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
+        Input* nux_input          = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
+        Input* nuy_input          = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
+        Input* potential_input    = element->GetInput(PotentialEnum);             _assert_(potential_input);
+        Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
 
 
@@ -326 +337 @@
                 nuy_input->GetInputValue(&nuy,gauss);
                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+                thicknessobs_input->GetInputValue(&thicknessobs,gauss);
+
+                dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
+                nuvobs2   = nux*vxobs*nux*vxobs + nuy*vyobs*nuy*vyobs;
+
 
                 /*Build gradient vector (actually -dJ/da): */
                 for(int i=0;i<numvertices;i++){
-                        ge[i]+= - weight*Jdet*gauss->weight*(+vyobs)*dphi[0]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
+                        //ge[i]+= - weight*Jdet*gauss->weight*(+vyobs)*dphi[0]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
+                        ge[i]+= - weight*Jdet*gauss->weight*basis[i]*(
+                                                -2.*nuy*vyobs*vyobs*(
+                                                        dphinorm2*dphinorm2 - thicknessobs*thicknessobs*dphinorm2*nuvobs2
+                                                        )/ ( (pow(nux,6)+pow(nuy,6))*pow(vyobs,6) + 3*(nux*vxobs*nuy*vyobs)*(nux*vxobs*nuy*vyobs)*nuvobs2)
+                                                );
                         _assert_(!xIsNan<IssmDouble>(ge[i]));
                 }

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

@@ -1355 +1355 @@
 
                         /*J = (H^2 - Hobs^2)^2*/
-                        //J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
+                        J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
                         /*J = phi^2*/
                         //J +=.5*potential*potential*weight*Jdet*gauss->weight;// OK
@@ -1363 +1363 @@
                         //J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
                         /*J = 1/2 (vbar ^ gard(phi))^2*/
-                        J += 0.5*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
+                        //J += 0.5*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
                 }