Changeset 11310


Ignore:
Timestamp:
02/02/12 13:59:41 (13 years ago)
Author:
Mathieu Morlighem
Message:

Simplifying gradient call in CM

Location:
issm/trunk-jpl/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp

    r9719 r11310  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void Gradjx( Vec* pgradient, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int control_type){
     12void Gradjx( Vec* pgradient,double* pnorm_grad, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int control_index){
    1313
    14         int  i;
    15         int  dim;
    16         int  numberofvertices;
    17         Vec  gradient = NULL;
     14        int    i;
     15        int    numberofvertices;
     16        double norm_grad;
     17        int   *control_type   = NULL;
     18        Vec    gradient = NULL;
    1819       
    1920        /*retrieve some parameters: */
    20         parameters->FindParam(&dim,MeshDimensionEnum);
     21        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    2122        numberofvertices=vertices->NumberOfVertices();
    2223
    … …  
    2728        for (i=0;i<elements->Size();i++){
    2829                Element* element=(Element*)elements->GetObjectByOffset(i);
    29                 element->Gradj(gradient,control_type);
     30                element->Gradj(gradient,control_type[control_index]);
    3031        }
    3132
    … …  
    3435        VecAssemblyEnd(gradient);
    3536
    36         /*Assign output pointers: */
     37        /*Clean up and return*/
     38        xfree((void**)&control_type);
     39        if(pnorm_grad){
     40                VecNorm(gradient,NORM_INFINITY,&norm_grad);
     41                *pnorm_grad=norm_grad;
     42        }
    3743        *pgradient=gradient;
    3844}

  • issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h

    r4236 r11310  
    1010
    1111/* local prototypes: */
    12 void Gradjx(Vec* pgrad_g, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,  Parameters* parameters, int control_type);
     12void Gradjx(Vec* pgrad_g,double* pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,  Parameters* parameters, int control_type);
    1313
    1414#endif  /* _GRADJX_H */

  • issm/trunk-jpl/src/c/solutions/controltao_core.cpp

    r11297 r11310  
    8585        adjointdiagnostic_core(user->femmodel);
    8686        //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum);
    87         Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
     87        Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,0);
    8888        VecScale(gradient,-1.);
    8989        VecCopy(gradient,G);

  • issm/trunk-jpl/src/c/solutions/gradient_core.cpp

    r11306 r11310  
    2020        int    *control_type   = NULL;
    2121        double *optscal_list   = NULL;
    22         double  optscal,norm_grad;
     22        double  optscal=1.e+100,norm_grad;
    2323
    2424        /*Intermediaries*/
    … …  
    3737
    3838                _printf_(VerboseControl(),"   compute gradient of J with respect to %s\n",EnumToStringx(control_type[i]));
    39                 Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type[i]);
     39                Gradjx(&gradient,&norm_grad,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,i);
    4040
    4141                if (orthogonalize){
    … …  
    5252                if(norm_grad<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm of J with respect to %s is zero",EnumToStringx(control_type[i]));
    5353                if(isnan(norm_grad))_error_("||∂J/∂α||∞ = NaN  gradient norm of J with respect to %s is NaN" ,EnumToStringx(control_type[i]));
    54                 optscal=optscal_list[num_controls*step+i]/norm_grad;
     54                optscal=min(optscal,optscal_list[num_controls*step+i]/norm_grad);
    5555
    5656                /*plug back into inputs: */
    57                 ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,control_type[i],new_gradient);
     57                ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],new_gradient);
    5858                VecFree(&new_gradient);
    5959        }
    6060
    6161        /*Scale Gradients*/
    62         for (int i=0;i<num_controls;i++) ControlInputScaleGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,control_type[i],optscal);
     62        for (int i=0;i<num_controls;i++) ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],optscal);
    6363
    6464        /*Clean up and return*/

  • issm/trunk-jpl/src/m/solutions/gradient_core.m

    r11306 r11310  
    2020end
    2121
    22         %recover parameters common to all solutions
    23         num_controls=femmodel.parameters.InversionNumControlParameters;
    24         control_type=femmodel.parameters.InversionControlParameters;
    25         control_steady=femmodel.parameters.ControlSteady;
    26         gradient_scaling_list=femmodel.parameters.InversionGradientScaling;
     22%recover parameters common to all solutions
     23num_controls=femmodel.parameters.InversionNumControlParameters;
     24control_type=femmodel.parameters.InversionControlParameters;
     25control_steady=femmodel.parameters.ControlSteady;
     26gradient_scaling_list=femmodel.parameters.InversionGradientScaling;
     27gradient_scaling=10^100;
    2728
    28         for i=1:num_controls,
     29for i=1:num_controls,
    2930
    30                 issmprintf(VerboseControl,['   compute gradient of J with respect to %s'],EnumToString(control_type(i)));
    31                 grad=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type(i));
     31        issmprintf(VerboseControl,['   compute gradient of J with respect to %s'],EnumToString(control_type(i)));
     32        [grad norm_grad]=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,i-1);
    3233
    33                 if orthogonalize,
    34                         issmprintf(VerboseControl,'%s',['   orthogonalization']);
    35                         old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,control_type(i));
    36                         new_gradient=Orth(grad,old_gradient);
    37                 else
    38                         new_gradient=grad;
    39                 end
    40 
    41                  %Get scaling factor of current control:
    42                  norm_grad=norm(new_gradient,inf);
    43                  if(norm_grad<=0),     error(['||∂J/∂α||∞ = 0   gradient norm of J with respect to ' EnumToString(control_type(i))  ' is zero']); end
    44                  if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i))  ' is NaN' ]); end
    45                  gradient_scaling=gradient_scaling_list(step,i)/norm_grad;
    46 
    47                 %plug back into inputs:
    48                 [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);
     34        if orthogonalize,
     35                issmprintf(VerboseControl,'%s',['   orthogonalization']);
     36                old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,control_type(i));
     37                new_gradient=Orth(grad,old_gradient);
     38        else
     39                new_gradient=grad;
    4940        end
    5041
    51         %Scale all gradients
    52         for i=1:num_controls,
    53                 [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);
    54         end
     42         %Get scaling factor of current control:
     43         norm_grad=norm(new_gradient,inf);
     44         if(norm_grad<=0),     error(['||∂J/∂α||∞ = 0   gradient norm of J with respect to ' EnumToString(control_type(i))  ' is zero']); end
     45         if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i))  ' is NaN' ]); end
     46         gradient_scaling=min(gradient_scaling,gradient_scaling_list(step,i)/norm_grad);
     47
     48        %plug back into inputs:
     49        [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);
     50
     51end
     52
     53%Scale all gradients
     54for i=1:num_controls,
     55        [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);
     56end

  • issm/trunk-jpl/src/mex/Gradj/Gradj.cpp

    r8910 r11310  
    88
    99        /*input datasets: */
    10         int         control_type;
     10        int         control_index;
     11        double      grad_norm;
    1112        Elements   *elements     = NULL;
    1213        Nodes      *nodes        = NULL;
    … …  
    3233        FetchMatlabData((DataSet**)&materials,MATERIALS);
    3334        FetchMatlabData(&parameters,PARAMETERS);
    34         FetchMatlabData(&control_type,CONTROLTYPE);
     35        FetchMatlabData(&control_index,CONTROLINDEX);
    3536
    3637        /*configure: */
    … …  
    4041
    4142        /*!Call core code: */
    42         Gradjx(&gradient, elements,nodes, vertices,loads, materials,parameters, control_type);
     43        Gradjx(&gradient,&grad_norm,elements,nodes, vertices,loads, materials,parameters, control_index);
    4344
    4445        /*write output : */
     46        WriteMatlabData(GRADNORM,grad_norm);
    4547        WriteMatlabData(GRADG,gradient);
    4648

  • issm/trunk-jpl/src/mex/Gradj/Gradj.h

    r6261 r11310  
    2424#define MATERIALS (mxArray*)prhs[4]
    2525#define PARAMETERS (mxArray*)prhs[5]
    26 #define CONTROLTYPE (mxArray*)prhs[6]
     26#define CONTROLINDEX (mxArray*)prhs[6]
    2727
    2828/* serial output macros: */
    2929#define GRADG (mxArray**)&plhs[0]
     30#define GRADNORM (mxArray**)&plhs[1]
    3031
    3132/* serial arg counts: */
    3233#undef NLHS
    33 #define NLHS  1
     34#define NLHS  2
    3435#undef NRHS
    3536#define NRHS  7
Note: See TracChangeset for help on using the changeset viewer.