Changeset 2355

Timestamp:

09/29/09 16:42:37 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added max and min dampening in misfit and gradient

File:

• issm/trunk/src/c/objects/Tria.cpp (modified) (10 diffs)

issm/trunk/src/c/objects/Tria.cpp

@@ -2394 +2394 @@
                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
                 for (i=0;i<numgrids;i++){
-                        grade_g_gaussian[i]=
-                          -2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i]                         //standard term dJ/dki
-                          -numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);  // regularization term d/dki(1/2*(dk/dx)^2)
+
+                        //standard term dJ/dki
+                        grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
+
+                        //noise dampening d/dki(1/2*(dk/dx)^2)
+                        grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+                        
+                        //min dampening
+                        if(drag<numpar->cm_mindmp_value){ 
+                                grade_g_gaussian[i]+=numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        }
+
+                        //max dampening
+                        if(drag<numpar->cm_maxdmp_value){ 
+                                grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        }
                 }
                 
@@ -2619 +2632 @@
 
                         //Add regularization term
-                        grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+                        grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+
+                        //min dampening
+                        if(drag<numpar->cm_mindmp_value){ 
+                                grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        }
+
+                        //max dampening
+                        if(drag<numpar->cm_maxdmp_value){ 
+                                grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        }
                 }
 
@@ -2689 +2712 @@
         double adjy_list[numgrids];
         double adjxadjy_list[numgrids][2];
+        double B[numgrids];
 
         int    dofs1[1]={0};
@@ -2725 +2749 @@
         double  thickness;
         int     dofs[1]={0};
-        double  dk[NDOF2]; 
+        double  dB[NDOF2]; 
+        double  B_gauss;
         
         ParameterInputs* inputs=NULL;
@@ -2747 +2772 @@
         if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
                 throw ErrorException(__FUNCT__,"missing adjoint input parameter");
+        }
+        if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){
+                throw ErrorException(__FUNCT__,"parameter B not found in input");
         }
 
@@ -2800 +2828 @@
                 GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
 
-                /*Get k derivative: dk/dx */
-                GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
+                /*Get B derivative: dB/dx */
+                GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3);
+                GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3);
 
                 /*Build gradje_g_gaussian vector (actually -dJ/dB): */
@@ -2809 +2838 @@
 
                         //Add regularization term
-                        grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+                        grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dB[0]+dh1dh2dh3_basic[1][i]*dB[1]);
+
+                        //min dampening
+                        if(B_gauss<numpar->cm_mindmp_value){ 
+                                grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        }
+
+                        //max dampening
+                        if(B_gauss<numpar->cm_maxdmp_value){ 
+                                grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        }
+
                 }
 
@@ -2866 +2906 @@
         double  gauss_weight;
         double  gauss_l1l2l3[3];
+        double  k_gauss;
+        double  B_gauss;
 
         /* parameters: */
@@ -2963 +3005 @@
                 #endif
                 
-                /*Add dampening term to misfit*/
+                /*Add dampening terms to misfit*/
                 if (strcmp(numpar->control_type,"drag")==0){
                         if (!shelf){
+
+                                //noise dampening
                                 GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3);
                                 Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
+
+                                //min dampening
+                                GetParameterValue(&k_gauss, &k[0],gauss_l1l2l3);
+                                if(k_gauss<numpar->cm_mindmp_value){ 
+                                        Jelem+=numpar->cm_mindmp_slope*k_gauss*Jdet*gauss_weight;
+                                }
+
+                                //max dampening
+                                if(k_gauss>numpar->cm_maxdmp_value){ 
+                                        Jelem+=numpar->cm_maxdmp_slope*k_gauss*Jdet*gauss_weight;
+                                }
                         }
                 }
@@ -2974 +3029 @@
                                 throw ErrorException(__FUNCT__,"parameter B not found in input");
                         }
+                        //noise dampening
                         GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3);
                         Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
+
+                        //min dampening
+                        GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3);
+                        if(B_gauss<numpar->cm_mindmp_value){ 
+                                Jelem+=numpar->cm_mindmp_slope*B_gauss*Jdet*gauss_weight;
+                        }
+
+                        //max dampening
+                        if(B_gauss>numpar->cm_maxdmp_value){ 
+                                Jelem+=numpar->cm_maxdmp_slope*B_gauss*Jdet*gauss_weight;
+                        }
                 }
                 else{