Index: /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 6225) +++ /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 6226) @@ -43,6 +43,5 @@ parameters->AddObject(new DoubleParam(MeanVelEnum,iomodel->meanvel)); parameters->AddObject(new DoubleParam(CmNoiseDmpEnum,iomodel->cm_noisedmp)); - parameters->AddObject(new DoubleParam(CmMinEnum,iomodel->cm_min)); - parameters->AddObject(new DoubleParam(CmMaxEnum,iomodel->cm_max)); + parameters->AddObject(new BoolParam(CmGradientEnum,iomodel->cm_gradient)); @@ -50,4 +49,6 @@ IoModelFetchData(&iomodel->cm_responses,NULL,NULL,iomodel_handle,"cm_responses"); IoModelFetchData(&iomodel->cm_jump,NULL,NULL,iomodel_handle,"cm_jump"); + IoModelFetchData(&iomodel->cm_min,NULL,NULL,iomodel_handle,"cm_min"); + IoModelFetchData(&iomodel->cm_max,NULL,NULL,iomodel_handle,"cm_max"); IoModelFetchData(&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal"); IoModelFetchData(&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter"); @@ -55,9 +56,13 @@ parameters->AddObject(new DoubleVecParam(CmResponsesEnum,iomodel->cm_responses,iomodel->nsteps)); parameters->AddObject(new DoubleVecParam(CmJumpEnum,iomodel->cm_jump,iomodel->nsteps)); - parameters->AddObject(new DoubleVecParam(OptScalEnum,iomodel->optscal,iomodel->nsteps)); + parameters->AddObject(new DoubleMatParam(OptScalEnum,iomodel->optscal,iomodel->nsteps,iomodel->num_control_type)); + parameters->AddObject(new DoubleVecParam(CmMinEnum,iomodel->cm_min,iomodel->num_control_type)); + parameters->AddObject(new DoubleVecParam(CmMaxEnum,iomodel->cm_max,iomodel->num_control_type)); parameters->AddObject(new DoubleVecParam(MaxIterEnum,iomodel->maxiter,iomodel->nsteps)); xfree((void**)&iomodel->cm_responses); xfree((void**)&iomodel->cm_jump); + xfree((void**)&iomodel->cm_min); + xfree((void**)&iomodel->cm_max); xfree((void**)&iomodel->optscal); xfree((void**)&iomodel->maxiter); Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6225) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6226) @@ -1134,9 +1134,10 @@ int* control_type=NULL; Input* input=NULL; - double cm_min,cm_max; + double *cm_min=NULL; + double *cm_max=NULL; /*retrieve some parameters: */ - this->parameters->FindParam(&cm_min,CmMinEnum); - this->parameters->FindParam(&cm_max,CmMaxEnum); + this->parameters->FindParam(&cm_min,NULL,CmMinEnum); + this->parameters->FindParam(&cm_max,NULL,CmMaxEnum); this->parameters->FindParam(&num_controls,NumControlsEnum); this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); @@ -1155,5 +1156,5 @@ ((ControlInput*)input)->UpdateValue(scalar); - input->Constrain(cm_min,cm_max); + input->Constrain(cm_min[i],cm_max[i]); if (save_parameter) ((ControlInput*)input)->SaveValue(); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6225) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6226) @@ -1342,9 +1342,10 @@ int* control_type=NULL; Input* input=NULL; - double cm_min,cm_max; + double *cm_min=NULL; + double *cm_max=NULL; /*retrieve some parameters: */ - this->parameters->FindParam(&cm_min,CmMinEnum); - this->parameters->FindParam(&cm_max,CmMaxEnum); + this->parameters->FindParam(&cm_min,NULL,CmMinEnum); + this->parameters->FindParam(&cm_max,NULL,CmMaxEnum); this->parameters->FindParam(&num_controls,NumControlsEnum); this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); @@ -1364,5 +1365,5 @@ ((ControlInput*)input)->UpdateValue(scalar); - input->Constrain(cm_min,cm_max); + input->Constrain(cm_min[i],cm_max[i]); if (save_parameter) ((ControlInput*)input)->SaveValue(); Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 6225) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 6226) @@ -88,4 +88,6 @@ xfree((void**)&this->weights); xfree((void**)&this->cm_jump); + xfree((void**)&this->cm_min); + xfree((void**)&this->cm_max); xfree((void**)&this->optscal); xfree((void**)&this->maxiter); @@ -166,6 +168,4 @@ IoModelFetchData(&this->tolx,iomodel_handle,"tolx"); IoModelFetchData(&this->cm_noisedmp,iomodel_handle,"cm_noisedmp"); - IoModelFetchData(&this->cm_min,iomodel_handle,"cm_min"); - IoModelFetchData(&this->cm_max,iomodel_handle,"cm_max"); IoModelFetchData(&this->cm_gradient,iomodel_handle,"cm_gradient"); IoModelFetchData(&this->eps_res,iomodel_handle,"eps_res"); @@ -306,6 +306,6 @@ this->maxiter=NULL; this->cm_noisedmp=0; - this->cm_min=0; - this->cm_max=0; + this->cm_min=NULL; + this->cm_max=NULL; this->cm_gradient=0; this->verbose=0; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 6225) +++ /issm/trunk/src/c/objects/IoModel.h (revision 6226) @@ -114,6 +114,6 @@ double stokesreconditioning; double cm_noisedmp; - double cm_min; - double cm_max; + double* cm_min; + double* cm_max; int cm_gradient;; Index: /issm/trunk/src/c/objects/Params/DoubleParam.cpp =================================================================== --- /issm/trunk/src/c/objects/Params/DoubleParam.cpp (revision 6225) +++ /issm/trunk/src/c/objects/Params/DoubleParam.cpp (revision 6226) @@ -167,4 +167,37 @@ } /*}}}*/ +/*FUNCTION DoubleParam::GetParameterValue(double** pdoublearray,int* pM){{{1*/ +void DoubleParam::GetParameterValue(double** pdoublearray,int* pM){ +#ifdef _SERIAL_ + double* output=NULL; + + output=(double*)xmalloc(1*sizeof(double)); + *output=(double)value; + + /*Assign output podoubleers:*/ + if(pM) *pM=1; + *pdoublearray=output; +#else + ISSMERROR("Double param of enum %i (%s) cannot return an array of double",enum_type,EnumToString(enum_type)); +#endif +} +/*}}}*/ +/*FUNCTION DoubleParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){{{1*/ +void DoubleParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){ +#ifdef _SERIAL_ + double* output=NULL; + + output=(double*)xmalloc(1*sizeof(double)); + *output=(double)value; + + /*Assign output podoubleers:*/ + if(pM) *pM=1; + if(pN) *pN=1; + *pdoublearray=output; +#else + ISSMERROR("Double param of enum %i (%s) cannot return an array of double",enum_type,EnumToString(enum_type)); +#endif +} +/*}}}*/ /*FUNCTION DoubleParam::SetMatlabField{{{1*/ #ifdef _SERIAL_ Index: /issm/trunk/src/c/objects/Params/DoubleParam.h =================================================================== --- /issm/trunk/src/c/objects/Params/DoubleParam.h (revision 6225) +++ /issm/trunk/src/c/objects/Params/DoubleParam.h (revision 6226) @@ -56,6 +56,6 @@ void GetParameterValue(char** pstring){ISSMERROR("Double param of enum %i (%s) cannot return a string",enum_type,EnumToString(enum_type));} void GetParameterValue(char*** pstringarray,int* pM){ISSMERROR("Double param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToString(enum_type));} - void GetParameterValue(double** pdoublearray,int* pM){ISSMERROR("Double param of enum %i (%s) cannot return a double array",enum_type,EnumToString(enum_type));} - void GetParameterValue(double** pdoublearray,int* pM, int* pN){ISSMERROR("Double param of enum %i (%s) cannot return a double array",enum_type,EnumToString(enum_type));} + void GetParameterValue(double** pdoublearray,int* pM); + void GetParameterValue(double** pdoublearray,int* pM, int* pN); void GetParameterValue(double*** parray, int* pM,int** pmdims, int** pndims){ISSMERROR("Double param of enum %i (%s) cannot return a matrix array",enum_type,EnumToString(enum_type));} void GetParameterValue(Vec* pvec){ISSMERROR("Double param of enum %i (%s) cannot return a Vec",enum_type,EnumToString(enum_type));} Index: /issm/trunk/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/control_core.cpp (revision 6225) +++ /issm/trunk/src/c/solutions/control_core.cpp (revision 6226) @@ -30,5 +30,4 @@ int* control_type = NULL; double* responses=NULL; - double* optscal=NULL; double* maxiter=NULL; double* cm_jump=NULL; @@ -52,5 +51,4 @@ femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum); femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); - femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum); femmodel->parameters->FindParam(&cm_jump,NULL,CmJumpEnum); @@ -107,5 +105,5 @@ _printf_("%s\n"," updating parameter using optimized search scalar..."); //true means update save controls - InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true); + InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true); if(controlconvergence(J,responses,eps_cm,n)) break; @@ -133,5 +131,4 @@ xfree((void**)&control_type); xfree((void**)&responses); - xfree((void**)&optscal); xfree((void**)&maxiter); xfree((void**)&cm_jump); Index: /issm/trunk/src/c/solutions/gradient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/gradient_core.cpp (revision 6225) +++ /issm/trunk/src/c/solutions/gradient_core.cpp (revision 6226) @@ -16,8 +16,9 @@ /*parameters: */ - bool control_steady; - int num_controls; - int* control_type=NULL; - int verbose; + bool control_steady; + int verbose; + int num_controls; + int *control_type = NULL; + double *optscal_list = NULL; /*Intermediaries*/ @@ -30,4 +31,5 @@ femmodel->parameters->FindParam(&num_controls,NumControlsEnum); femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum); + femmodel->parameters->FindParam(&optscal_list,NULL,NULL,OptScalEnum); femmodel->parameters->FindParam(&verbose,VerboseEnum); @@ -51,4 +53,7 @@ VecFree(&gradient); + /*Scale gradient for current step and current parameter*/ + VecScale(new_gradient,optscal_list[num_controls*step+i]); + /*plug back into inputs: */ ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,control_type[i],new_gradient); Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 6225) +++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 6226) @@ -30,5 +30,4 @@ FemModel *femmodel = NULL; int n; - double *optscal = NULL; double *responses = NULL; int solution_type; @@ -43,5 +42,4 @@ n=optargs->n; - femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); femmodel->parameters->FindParam(&isstokes,IsStokesEnum); @@ -61,5 +59,5 @@ /*update parameter according to scalar: */ //false means: do not save control - InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],false); + InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false); /*Run diagnostic with updated inputs: */ @@ -82,5 +80,4 @@ /*Free ressources:*/ xfree((void**)&responses); - xfree((void**)&optscal); return J; } Index: /issm/trunk/src/m/model/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 6225) +++ /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 6226) @@ -170,9 +170,4 @@ end %}}} -%SCALAR {{{ -if ~isscalar(md.control_type), - error('model not consistent: md.control_type should be a scalar'); -end -%}}} %PARAMETEROUTPUT {{{1 if md.numoutput~=length(md.parameteroutput), @@ -344,12 +339,14 @@ %CONTROL TYPE - if ~isnumeric(md.control_type), - error('model not consistent: control_type should be an enum'); - end + num_controls=numel(md.control_type); checkvalues(md,{'control_type'},[DhDtEnum DragCoefficientEnum RheologyBbarEnum VxEnum VyEnum]); %LENGTH CONTROL FIELDS - fields={'maxiter','optscal','cm_responses','cm_jump'}; + fields={'maxiter','cm_responses','cm_jump'}; checksize(md,fields,[md.nsteps 1]); + fields={'optscal'}; + checksize(md,fields,[md.nsteps num_controls]); + fields={'cm_min','cm_max'}; + checksize(md,fields,[num_controls 1]); %RESPONSES Index: /issm/trunk/src/m/model/marshall.m =================================================================== --- /issm/trunk/src/m/model/marshall.m (revision 6225) +++ /issm/trunk/src/m/model/marshall.m (revision 6226) @@ -119,6 +119,6 @@ WriteData(fid,md.maxiter,'Mat','maxiter'); WriteData(fid,md.cm_noisedmp,'Scalar','cm_noisedmp'); -WriteData(fid,md.cm_min,'Scalar','cm_min'); -WriteData(fid,md.cm_max,'Scalar','cm_max'); +WriteData(fid,md.cm_min,'Mat','cm_min'); +WriteData(fid,md.cm_max,'Mat','cm_max'); WriteData(fid,md.cm_gradient,'Integer','cm_gradient'); WriteData(fid,md.eps_res,'Scalar','eps_res'); Index: /issm/trunk/src/m/solutions/control_core.m =================================================================== --- /issm/trunk/src/m/solutions/control_core.m (revision 6225) +++ /issm/trunk/src/m/solutions/control_core.m (revision 6226) @@ -13,5 +13,4 @@ solution_type=femmodel.parameters.SolutionType; responses=femmodel.parameters.CmResponses; - optscal=femmodel.parameters.OptScal; maxiter=femmodel.parameters.MaxIter; cm_jump=femmodel.parameters.CmJump; @@ -66,5 +65,5 @@ displaystring('\n%s',[' updating parameter using optimized search scalar:']); - [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal(n),1); + [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar,1); disp([' value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]); Index: /issm/trunk/src/m/solutions/gradient_core.m =================================================================== --- /issm/trunk/src/m/solutions/gradient_core.m (revision 6225) +++ /issm/trunk/src/m/solutions/gradient_core.m (revision 6226) @@ -25,4 +25,5 @@ control_type=femmodel.parameters.ControlType; control_steady=femmodel.parameters.ControlSteady; + optscal=femmodel.parameters.OptScal; for i=1:num_controls, @@ -45,4 +46,7 @@ end + %Scale Gradient + new_gradient=optscal(step,i)*new_gradient; + %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); Index: /issm/trunk/src/m/solutions/objectivefunctionC.m =================================================================== --- /issm/trunk/src/m/solutions/objectivefunctionC.m (revision 6225) +++ /issm/trunk/src/m/solutions/objectivefunctionC.m (revision 6226) @@ -4,5 +4,4 @@ conserve_loads=true; %recover some parameters -optscal=femmodel.parameters.OptScal(n); response=femmodel.parameters.CmResponses(n); analysis_type=femmodel.parameters.AnalysisType; @@ -20,5 +19,5 @@ %Use search scalar to shoot parameter in the gradient direction: -[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal,0); +[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar,0); %Run diagnostic with updated inputs: