Changeset 17905

Timestamp:

05/01/14 10:17:42 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added constraints on m1qn3 inversion

File:

• issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp

@@ -35 +35 @@
         IssmDouble  *X  = NULL;
         IssmDouble  *G  = NULL;
-        IssmDouble  *XL = NULL;
-        IssmDouble  *XU = NULL;
 
         /*Recover some parameters*/
@@ -139 +137 @@
         femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
 
+        /*Constrain input vector*/
+        IssmDouble  *XL = NULL;
+        IssmDouble  *XU = NULL;
+        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+        for(long i=0;i<*n;i++){
+                if(X[i]>XU[i]) X[i]=XU[i];
+                if(X[i]<XL[i]) X[i]=XL[i];
+        }
+
         /*Update control input*/
         SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
@@ -162 +170 @@
         _printf0_("\n");
 
-        if(indic==0) return; /*dry run, no gradient required*/
+        if(indic==0){
+                /*dry run, no gradient required*/
+                xDelete<int>(responses);
+                xDelete<IssmDouble>(XU);
+                xDelete<IssmDouble>(XL);
+                return;
+        }
 
         /*Compute Adjoint*/
@@ -172 +186 @@
         Gradjx(&G2,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
         for(long i=0;i<*n;i++) G[i] = -G2[i];
+        xDelete<IssmDouble>(G2);
+
+        /*Constrain Gradient*/
+        for(long i=0;i<*n;i++){
+                if(X[i]>=XU[i]) G[i]=0.;
+                if(X[i]<=XL[i]) G[i]=0.;
+        }
 
         /*Clean-up and return*/
         xDelete<int>(responses);
-        xDelete<IssmDouble>(G2);
+        xDelete<IssmDouble>(XU);
+        xDelete<IssmDouble>(XL);
 }