Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 11309) +++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 11310) @@ -10,13 +10,14 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void Gradjx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int control_type){ +void Gradjx( Vec* pgradient,double* pnorm_grad, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int control_index){ - int i; - int dim; - int numberofvertices; - Vec gradient = NULL; + int i; + int numberofvertices; + double norm_grad; + int *control_type = NULL; + Vec gradient = NULL; /*retrieve some parameters: */ - parameters->FindParam(&dim,MeshDimensionEnum); + parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); numberofvertices=vertices->NumberOfVertices(); @@ -27,5 +28,5 @@ for (i=0;iSize();i++){ Element* element=(Element*)elements->GetObjectByOffset(i); - element->Gradj(gradient,control_type); + element->Gradj(gradient,control_type[control_index]); } @@ -34,5 +35,10 @@ VecAssemblyEnd(gradient); - /*Assign output pointers: */ + /*Clean up and return*/ + xfree((void**)&control_type); + if(pnorm_grad){ + VecNorm(gradient,NORM_INFINITY,&norm_grad); + *pnorm_grad=norm_grad; + } *pgradient=gradient; } Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 11309) +++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 11310) @@ -10,5 +10,5 @@ /* local prototypes: */ -void Gradjx(Vec* pgrad_g, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int control_type); +void Gradjx(Vec* pgrad_g,double* pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int control_type); #endif /* _GRADJX_H */ Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 11309) +++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp (revision 11310) @@ -85,5 +85,5 @@ adjointdiagnostic_core(user->femmodel); //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum); - Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum); + Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,0); VecScale(gradient,-1.); VecCopy(gradient,G); Index: /issm/trunk-jpl/src/c/solutions/gradient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/gradient_core.cpp (revision 11309) +++ /issm/trunk-jpl/src/c/solutions/gradient_core.cpp (revision 11310) @@ -20,5 +20,5 @@ int *control_type = NULL; double *optscal_list = NULL; - double optscal,norm_grad; + double optscal=1.e+100,norm_grad; /*Intermediaries*/ @@ -37,5 +37,5 @@ _printf_(VerboseControl()," compute gradient of J with respect to %s\n",EnumToStringx(control_type[i])); - Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type[i]); + Gradjx(&gradient,&norm_grad,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,i); if (orthogonalize){ @@ -52,13 +52,13 @@ if(norm_grad<=0) _error_("||∂J/∂α||∞ = 0 gradient norm of J with respect to %s is zero",EnumToStringx(control_type[i])); if(isnan(norm_grad))_error_("||∂J/∂α||∞ = NaN gradient norm of J with respect to %s is NaN" ,EnumToStringx(control_type[i])); - optscal=optscal_list[num_controls*step+i]/norm_grad; + optscal=min(optscal,optscal_list[num_controls*step+i]/norm_grad); /*plug back into inputs: */ - ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,control_type[i],new_gradient); + ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],new_gradient); VecFree(&new_gradient); } /*Scale Gradients*/ - for (int i=0;i elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,control_type[i],optscal); + for (int i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],optscal); /*Clean up and return*/ Index: /issm/trunk-jpl/src/m/solutions/gradient_core.m =================================================================== --- /issm/trunk-jpl/src/m/solutions/gradient_core.m (revision 11309) +++ /issm/trunk-jpl/src/m/solutions/gradient_core.m (revision 11310) @@ -20,35 +20,37 @@ end - %recover parameters common to all solutions - num_controls=femmodel.parameters.InversionNumControlParameters; - control_type=femmodel.parameters.InversionControlParameters; - control_steady=femmodel.parameters.ControlSteady; - gradient_scaling_list=femmodel.parameters.InversionGradientScaling; +%recover parameters common to all solutions +num_controls=femmodel.parameters.InversionNumControlParameters; +control_type=femmodel.parameters.InversionControlParameters; +control_steady=femmodel.parameters.ControlSteady; +gradient_scaling_list=femmodel.parameters.InversionGradientScaling; +gradient_scaling=10^100; - for i=1:num_controls, +for i=1:num_controls, - issmprintf(VerboseControl,[' compute gradient of J with respect to %s'],EnumToString(control_type(i))); - grad=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type(i)); + issmprintf(VerboseControl,[' compute gradient of J with respect to %s'],EnumToString(control_type(i))); + [grad norm_grad]=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,i-1); - if orthogonalize, - issmprintf(VerboseControl,'%s',[' orthogonalization']); - old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,control_type(i)); - new_gradient=Orth(grad,old_gradient); - else - new_gradient=grad; - end - - %Get scaling factor of current control: - norm_grad=norm(new_gradient,inf); - if(norm_grad<=0), error(['||∂J/∂α||∞ = 0 gradient norm of J with respect to ' EnumToString(control_type(i)) ' is zero']); end - if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i)) ' is NaN' ]); end - gradient_scaling=gradient_scaling_list(step,i)/norm_grad; - - %plug back into inputs: - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputSetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,control_type(i),new_gradient); + if orthogonalize, + issmprintf(VerboseControl,'%s',[' orthogonalization']); + old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,control_type(i)); + new_gradient=Orth(grad,old_gradient); + else + new_gradient=grad; end - %Scale all gradients - for i=1:num_controls, - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputScaleGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,control_type(i),gradient_scaling); - end + %Get scaling factor of current control: + norm_grad=norm(new_gradient,inf); + if(norm_grad<=0), error(['||∂J/∂α||∞ = 0 gradient norm of J with respect to ' EnumToString(control_type(i)) ' is zero']); end + if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i)) ' is NaN' ]); end + gradient_scaling=min(gradient_scaling,gradient_scaling_list(step,i)/norm_grad); + + %plug back into inputs: + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputSetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,control_type(i),new_gradient); + +end + +%Scale all gradients +for i=1:num_controls, + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputScaleGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials, femmodel.parameters,control_type(i),gradient_scaling); +end Index: /issm/trunk-jpl/src/mex/Gradj/Gradj.cpp =================================================================== --- /issm/trunk-jpl/src/mex/Gradj/Gradj.cpp (revision 11309) +++ /issm/trunk-jpl/src/mex/Gradj/Gradj.cpp (revision 11310) @@ -8,5 +8,6 @@ /*input datasets: */ - int control_type; + int control_index; + double grad_norm; Elements *elements = NULL; Nodes *nodes = NULL; @@ -32,5 +33,5 @@ FetchMatlabData((DataSet**)&materials,MATERIALS); FetchMatlabData(¶meters,PARAMETERS); - FetchMatlabData(&control_type,CONTROLTYPE); + FetchMatlabData(&control_index,CONTROLINDEX); /*configure: */ @@ -40,7 +41,8 @@ /*!Call core code: */ - Gradjx(&gradient, elements,nodes, vertices,loads, materials,parameters, control_type); + Gradjx(&gradient,&grad_norm,elements,nodes, vertices,loads, materials,parameters, control_index); /*write output : */ + WriteMatlabData(GRADNORM,grad_norm); WriteMatlabData(GRADG,gradient); Index: /issm/trunk-jpl/src/mex/Gradj/Gradj.h =================================================================== --- /issm/trunk-jpl/src/mex/Gradj/Gradj.h (revision 11309) +++ /issm/trunk-jpl/src/mex/Gradj/Gradj.h (revision 11310) @@ -24,12 +24,13 @@ #define MATERIALS (mxArray*)prhs[4] #define PARAMETERS (mxArray*)prhs[5] -#define CONTROLTYPE (mxArray*)prhs[6] +#define CONTROLINDEX (mxArray*)prhs[6] /* serial output macros: */ #define GRADG (mxArray**)&plhs[0] +#define GRADNORM (mxArray**)&plhs[1] /* serial arg counts: */ #undef NLHS -#define NLHS 1 +#define NLHS 2 #undef NRHS #define NRHS 7