Changeset 27724

Timestamp:

05/04/23 11:29:48 (23 months ago)

Author:

Mathieu Morlighem

Message:

CHG: making gradient of param output

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Params/ControlParam.cpp (modified) (3 diffs)
• classes/Params/Parameters.cpp (modified) (1 diff)
• classes/Params/Parameters.h (modified) (1 diff)
• cores/controladm1qn3_core.cpp (modified) (3 diffs)
• cores/transient_core.cpp (modified) (2 diffs)
• modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp (modified) (1 diff)
• modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -3498 +3498 @@
                 if(!IsOnBase()) return;
         }
-
-        /*Get out if this is not an element input*/
-        if(!IsInputEnum(control_enum)) return;
 
         /*Prepare index list*/

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -5032 +5032 @@
                 }
         }
-
-        /*Get out if this is not an element input*/
-        if(!IsInputEnum(control_enum)) return;
 
         /*Get list of ids for this element and this control*/

issm/trunk-jpl/src/c/classes/Params/ControlParam.cpp

@@ -62 +62 @@
         
                 savedvalue=NULL;
-                gradient=NULL;
+                gradient=xNewZeroInit<IssmDouble>(M);
         }
         else{
@@ -78 +78 @@
         xDelete<IssmDouble>(minvalue);
         xDelete<IssmDouble>(maxvalue);
+        xDelete<IssmDouble>(gradient);
+        xDelete<IssmDouble>(savedvalue);
         return;
 }
@@ -201 +203 @@
 void  ControlParam::SetGradient(IssmDouble* poutput, int in_M){/*{{{*/
 
+        _assert_(in_M==this->M);
+
         /*avoid leak: */
         xDelete<IssmDouble>(this->gradient);

issm/trunk-jpl/src/c/classes/Params/Parameters.cpp

@@ -649 +649 @@
 }
 /*}}}*/
-void   Parameters::SetControlParam(IssmDouble* IssmDoublearray,int M, int enum_type){/*{{{*/
-
-        Param* param=NULL;
-
-        /*first, figure out if the param has already been created: */
-        param=xDynamicCast<Param*>(this->FindParamObject(enum_type));
-
-        if(param) param->SetValue(IssmDoublearray,M); //already exists, just set it.
+void   Parameters::SetControlFromVector(IssmDouble* vector, int enum_type, int M, int N, int offset){/*{{{*/
+
+        /*first, figure out if the param has already been created: */
+        Param* param=NULL;
+        param=xDynamicCast<Param*>(this->FindParamObject(enum_type));
+
+        if(param) param->SetValue(&vector[offset], M);
+        else _error_("Param "<< EnumToStringx(enum_type) << " cannot setValue");
+}
+/*}}}*/
+void   Parameters::SetGradientFromVector(IssmDouble* vector, int enum_type, int M, int N, int offset){/*{{{*/
+
+        /*first, figure out if the param has already been created: */
+        Param* param=NULL;
+        param=xDynamicCast<Param*>(this->FindParamObject(enum_type));
+
+        if(param) param->SetGradient(&vector[offset], M);
         else _error_("Param "<< EnumToStringx(enum_type) << " cannot setValue");
 }

issm/trunk-jpl/src/c/classes/Params/Parameters.h

@@ -74 +74 @@
                 void  SetParam(FILE* fid,int enum_type);
                 void  SetParam(DataSet* dataset,int enum_type);
-                void  SetControlParam(IssmDouble* IssmDoublearray,int M, int enum_type);
+                void  SetControlFromVector(IssmDouble* array, int enum_type, int M, int N, int offset);
+                void  SetGradientFromVector(IssmDouble* array, int enum_type, int M, int N, int offset);
                 void  ControlParamSetGradient(IssmDouble* IssmDoublearray, int M, int enum_type);
                 void  GetVectorFromControl(Vector<IssmDouble>* vector,int control_enum,int control_index,int N,const char* data,int offset);

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

@@ -251 +251 @@
         #ifdef _HAVE_CODIPACK_
         if(checkpoint_frequency && solution_type == TransientSolutionEnum){
-                if(IsParamEnum(control_enum[0])){
-#ifdef _HAVE_AD_
-                        IssmDouble* aX=xNew<IssmDouble>(intn,"t");
-#else
-                        IssmDouble* aX=xNew<IssmDouble>(intn);
-#endif
-                        for(int i=0;i<intn;i++) {
-                                aX[i]=X[i];
-                        }
-                        femmodel->parameters->SetControlParam(aX,M[0],control_enum[0]); 
-                }
-                else{
-                        SetControlInputsFromVectorx(femmodel,X);
-                }
+                SetControlInputsFromVectorx(femmodel,X);
                 *pf = transient_ad(femmodel, G, &Jlist[(*Jlisti)*JlistN]);
         }
@@ -317 +304 @@
 
                 ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-                if(IsParamEnum(control_enum[0])){
-                        femmodel->parameters->SetControlParam(aX,M[0],control_enum[0]);
-                }
-                else{
-                        SetControlInputsFromVectorx(femmodel,aX);
-                }
+                SetControlInputsFromVectorx(femmodel,aX);
                 xDelete<IssmDouble>(aX);
 
@@ -673 +655 @@
         }
         
-        if(IsParamEnum(control_enum[0])){
-                //femmodel->parameters->ControlParamSetGradient(aG,M[0],control_enum[0]);
-                femmodel->parameters->SetControlParam(aX,M[0],control_enum[0]); 
-        }
-        else{
-                ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
-                SetControlInputsFromVectorx(femmodel,aX);
-        }
+        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
+        SetControlInputsFromVectorx(femmodel,aX);
         xDelete(aX);
 

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

@@ -392 +392 @@
         delete hdl_regin;
         if(my_rank==0) for(int i=0; i < Xsize; i++) tape_codi.registerInput(X[i]);
-        
-        if(IsParamEnum(control_enum[0])){
-                femmodel->parameters->SetControlParam(X,M[0],control_enum[0]);
-        }
-        else{
-                SetControlInputsFromVectorx(femmodel,X);
-        }
+        SetControlInputsFromVectorx(femmodel,X);
         
         IssmDouble J     = 0.;
@@ -461 +455 @@
                 /*Tell codipack that X is the independent*/
                 for(int i=0; i<Xsize; i++) tape_codi.registerInput(X[i]);
-                
-                if(IsParamEnum(control_enum[0])){
-                        femmodel->parameters->SetControlParam(X,M[0],control_enum[0]);
-                }
-                else{
-                        SetControlInputsFromVectorx(femmodel,X);
-                }
+                SetControlInputsFromVectorx(femmodel,X);
 
                 /*Get New state*/

issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp

@@ -25 +25 @@
         int offset = 0;
         for(int i=0;i<num_controls;i++){
-                for(Object* & object : elements->objects){
-                        Element* element=xDynamicCast<Element*>(object);
-                        element->ControlInputSetGradient(gradient,control_type[i],i,offset,M_all[i],N_all[i],interp_all[i]);
+                /*Is the control a Param?*/
+                if(IsParamEnum(control_type[i])){
+                        parameters->SetGradientFromVector(gradient, control_type[i], M_all[i], N_all[i], offset);
+                }
+                else if(IsInputEnum(control_type[i])){
+                        for(Object* & object : elements->objects){
+                                Element* element=xDynamicCast<Element*>(object);
+                                element->ControlInputSetGradient(gradient,control_type[i],i,offset,M_all[i],N_all[i],interp_all[i]);
+                        }
+                }
+                else{
+                        _error_("not supported yet");
                 }
                 offset+=M_all[i]*N_all[i];

issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp

@@ -22 +22 @@
         int offset = 0;
         for(int i=0;i<num_controls;i++){
-                for(Object* & object : femmodel->elements->objects){
-                        Element* element=xDynamicCast<Element*>(object);
-                        element->SetControlInputsFromVector(vector,control_type[i],i,offset,M[i],N[i]);
+                /*Is the control a Param?*/
+                if(IsParamEnum(control_type[i])){
+                        femmodel->parameters->SetControlFromVector(vector,control_type[i],M[i],N[i],offset);
+                }
+                else if(IsInputEnum(control_type[i])){
+                        for(Object* & object : femmodel->elements->objects){
+                                Element* element=xDynamicCast<Element*>(object);
+                                element->SetControlInputsFromVector(vector,control_type[i],i,offset,M[i],N[i]);
+                        }
+                }
+                else{
+                        _error_("not supported yet");
                 }
                 offset += M[i]*N[i];