Changeset 18369

Timestamp:

08/12/14 16:38:41 (11 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed H-Hobs cost function for balancethickness2

Location:

issm/trunk-jpl/src/c

Files:

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

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

@@ -100 +100 @@
                 hu2 = hobs*hobs*vbarobs2;
 
-                /*H^2 - Hobs^2*/
-                NUMxH2 = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
-                NUMyH2 = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
-                DENH2 = vbarobs2*vbarobs2+1.e-14;
-
-                /*Ubar-Ubar_obs*/
-                NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0];
-                NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1];
-                DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
-
                 /*Loop over all requested responses*/
                 for(int resp=0;resp<num_responses;resp++){
@@ -117 +107 @@
                                 case Balancethickness2MisfitEnum:
                                         /*J = (H^2 - Hobs^2)^2*/
-                                        //for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight; NOT WORKING
+                                        for(i=0;i<numnodes;i++){
+
+                                                /*H^2 - Hobs^2*/
+                                                NUMxH2 = 2.*dbasis[0*numnodes+i]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
+                                                NUMyH2 = 2.*dbasis[1*numnodes+i]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
+                                                DENH2 = vbarobs2*vbarobs2+1.e-14;
+
+                                                /*Ubar-Ubar_obs*/
+                                                NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0*numnodes+i];
+                                                NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1*numnodes+i];
+                                                DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
+
+                                                pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight;
+                                        }
                                         /*J = phi^2*/
-                                        //for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; OK
+                                        //for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; //OK
+                                        /*J = grad J ^2*/
+                                        //for(i=0;i<numnodes;i++) pe->values[i]+= (dphi[0]*dbasis[0*numnodes+i] + dphi[1]*dbasis[1*numnodes+i])*weight*Jdet*gauss->weight; //OK
                                         /*J = (ubar - nux*uobs)^2*/
-                                        for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
+                                        //for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
                                         break;
                                 default:

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

@@ -1304 +1304 @@
         IssmDouble J_sum;
 
-        IssmDouble  weight,thicknessobs,thickness,potential;
+        IssmDouble  weight,thicknessobs,thickness,potential,dpotential[2];
         IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nux,nuy;
         IssmDouble  Jdet;
@@ -1343 +1343 @@
                         thicknessobs_input->GetInputValue(&thicknessobs,gauss);
                         potential_input->GetInputValue(&potential,gauss);
+                        potential_input->GetInputDerivativeValue(&dpotential[0],xyz_list,gauss);
                         vxobs_input->GetInputValue(&vxobs,gauss);
                         vyobs_input->GetInputValue(&vyobs,gauss);
@@ -1354 +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
+                        //J +=.5*potential*potential*weight*Jdet*gauss->weight;// OK
+                        /*J = grad phi^2*/
+                        //J +=.5*(dpotential[0]*dpotential[0] + dpotential[1]*dpotential[1])*weight*Jdet*gauss->weight;
                         /*J = (ubar - nux*uobs)^2*/
-                        J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
+                        //J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
                 }