Changeset 9563
- Timestamp:
- 09/01/11 15:20:52 (14 years ago)
- Location:
- issm/trunk/src
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
r9558 r9563 456 456 TolxEnum, 457 457 OptscalEnum, 458 EpsvelEnum,459 MeanvelEnum,460 458 OutputfilenameEnum, 461 459 WaterfractionEnum, -
issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp
r9558 r9563 399 399 case TolxEnum : return "Tolx"; 400 400 case OptscalEnum : return "Optscal"; 401 case EpsvelEnum : return "Epsvel";402 case MeanvelEnum : return "Meanvel";403 401 case OutputfilenameEnum : return "Outputfilename"; 404 402 case WaterfractionEnum : return "Waterfraction"; -
issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
r9476 r9563 39 39 parameters->AddObject(iomodel->CopyConstantObject(NumCmResponsesEnum)); 40 40 parameters->AddObject(iomodel->CopyConstantObject(NstepsEnum)); 41 parameters->AddObject(iomodel->CopyConstantObject(TolxEnum));42 41 parameters->AddObject(iomodel->CopyConstantObject(EpsCmEnum)); 43 42 parameters->AddObject(iomodel->CopyConstantObject(CmGradientEnum)); 44 parameters->AddObject(iomodel->CopyConstantObject(MeanvelEnum));45 43 46 44 /*What solution type?*/ -
issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
r9515 r9563 35 35 parameters->AddObject(iomodel->CopyConstantObject(MaxNonlinearIterationsEnum)); 36 36 parameters->AddObject(iomodel->CopyConstantObject(MaxSteadystateIterationsEnum)); 37 parameters->AddObject(iomodel->CopyConstantObject(EpsvelEnum));38 37 parameters->AddObject(iomodel->CopyConstantObject(YtsEnum)); 39 38 parameters->AddObject(iomodel->CopyConstantObject(DtEnum)); -
issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp
r9558 r9563 397 397 else if (strcmp(name,"Tolx")==0) return TolxEnum; 398 398 else if (strcmp(name,"Optscal")==0) return OptscalEnum; 399 else if (strcmp(name,"Epsvel")==0) return EpsvelEnum;400 else if (strcmp(name,"Meanvel")==0) return MeanvelEnum;401 399 else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum; 402 400 else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r9451 r9563 1840 1840 int artdiff; 1841 1841 int i,j,ig,found=0; 1842 double Jdet,u,v,w,um,vm,wm ,epsvel;1842 double Jdet,u,v,w,um,vm,wm; 1843 1843 double gravity,rho_ice,rho_water; 1844 double epsvel=1.e-13; 1844 1845 double heatcapacity,thermalconductivity,dt; 1845 1846 double pressure,enthalpy; … … 1875 1876 this->parameters->FindParam(&dt,DtEnum); 1876 1877 this->parameters->FindParam(&artdiff,ArtificialDiffusivityEnum); 1877 this->parameters->FindParam(&epsvel,EpsvelEnum);1878 1878 Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); 1879 1879 Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); … … 2117 2117 int artdiff; 2118 2118 int i,j,ig,found=0; 2119 double Jdet,u,v,w,um,vm,wm,epsvel; 2119 double Jdet,u,v,w,um,vm,wm; 2120 double epsvel=1.e-13; 2120 2121 double gravity,rho_ice,rho_water; 2121 2122 double heatcapacity,thermalconductivity,dt; … … 2149 2150 this->parameters->FindParam(&dt,DtEnum); 2150 2151 this->parameters->FindParam(&artdiff,ArtificialDiffusivityEnum); 2151 this->parameters->FindParam(&epsvel,EpsvelEnum);2152 2152 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); 2153 2153 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); -
issm/trunk/src/c/objects/Elements/Tria.cpp
r9512 r9563 1510 1510 double Jdet; 1511 1511 double obs_velocity_mag,velocity_mag; 1512 double dux,duy,meanvel,epsvel; 1512 double dux,duy; 1513 double epsvel=1.e-13; 1514 double meanvel=1000/(365*24*3600); /*1000 m/yr*/ 1513 1515 double scalex=0,scaley=0,scale=0,S=0; 1514 1516 double vx,vy,vxobs,vyobs,weight; … … 1524 1526 this->parameters->FindParam(&num_responses,NumCmResponsesEnum); 1525 1527 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 1526 this->parameters->FindParam(&meanvel,MeanvelEnum);1527 this->parameters->FindParam(&epsvel,EpsvelEnum);1528 1528 Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); 1529 1529 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); … … 1686 1686 double Jdet; 1687 1687 double obs_velocity_mag,velocity_mag; 1688 double dux,duy,meanvel,epsvel; 1688 double dux,duy; 1689 double epsvel=1.e-13; 1690 double meanvel=1000/(365*24*3600); /*1000 m/yr*/ 1689 1691 double scalex=0,scaley=0,scale=0,S=0; 1690 1692 double vx,vy,vxobs,vyobs,weight; … … 1700 1702 this->parameters->FindParam(&num_responses,NumCmResponsesEnum); 1701 1703 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 1702 this->parameters->FindParam(&meanvel,MeanvelEnum);1703 this->parameters->FindParam(&epsvel,EpsvelEnum);1704 1704 Input* weights_input = inputs->GetInput(WeightsEnum); _assert_(weights_input); 1705 1705 Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); … … 4856 4856 int i,ig; 4857 4857 double Jelem=0; 4858 double meanvel, epsvel,misfit,Jdet; 4858 double misfit,Jdet; 4859 double epsvel=1.e-13; 4860 double meanvel=1000/(365*24*3600); /*1000 m/yr*/ 4859 4861 double velocity_mag,obs_velocity_mag; 4860 4862 double xyz_list[NUMVERTICES][3]; … … 4869 4871 4870 4872 /*Retrieve all inputs we will be needing: */ 4871 this->parameters->FindParam(&meanvel,MeanvelEnum);4872 this->parameters->FindParam(&epsvel,EpsvelEnum);4873 4873 Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); 4874 4874 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); … … 4922 4922 int fit=-1; 4923 4923 double Jelem=0, S=0; 4924 double meanvel, epsvel, misfit, Jdet; 4924 double epsvel=1.e-13; 4925 double meanvel=1000/(365*24*3600); /*1000 m/yr*/ 4926 double misfit, Jdet; 4925 4927 double vx,vy,vxobs,vyobs,weight; 4926 4928 double xyz_list[NUMVERTICES][3]; … … 4934 4936 4935 4937 /*Retrieve all inputs we will be needing: */ 4936 this->parameters->FindParam(&meanvel,MeanvelEnum);4937 this->parameters->FindParam(&epsvel,EpsvelEnum);4938 4938 Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); 4939 4939 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); … … 5011 5011 double Jelem=0; 5012 5012 double scalex=1,scaley=1; 5013 double meanvel, epsvel,misfit,Jdet; 5013 double misfit,Jdet; 5014 double epsvel=1.e-13; 5015 double meanvel=1000/(365*24*3600); /*1000 m/yr*/ 5014 5016 double vx,vy,vxobs,vyobs,weight; 5015 5017 double xyz_list[NUMVERTICES][3]; … … 5023 5025 5024 5026 /*Retrieve all inputs we will be needing: */ 5025 this->parameters->FindParam(&meanvel,MeanvelEnum);5026 this->parameters->FindParam(&epsvel,EpsvelEnum);5027 5027 Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); 5028 5028 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); -
issm/trunk/src/c/objects/OptPars.h
r2267 r9563 10 10 double xmin; 11 11 double xmax; 12 double tolerance;13 12 double cm_jump; 14 13 int maxiter; -
issm/trunk/src/c/shared/Numerics/BrentSearch.cpp
r9320 r9563 28 28 double xmax,xmin,xbest; 29 29 double x,x1,x2,xm; 30 double tol1,tol2,seps,tolerance; 30 double tol1,tol2,seps; 31 double tolerance=1.e-4; 31 32 int maxiter,iter; 32 33 bool loop=true,goldenflag; … … 35 36 xmin=optpars->xmin; 36 37 xmax=optpars->xmax; 37 tolerance=optpars->tolerance;38 38 maxiter=optpars->maxiter; 39 39 cm_jump=optpars->cm_jump; -
issm/trunk/src/c/shared/Numerics/OptimalSearch.cpp
r9320 r9563 24 24 25 25 /*tolerances: */ 26 double seps,tolerance; 26 double seps; 27 double tolerance=1.e-4; 27 28 int maxiter; 28 29 … … 35 36 x1 =optpars->xmin; 36 37 x2 =optpars->xmax; 37 tolerance=optpars->tolerance;38 38 maxiter =optpars->maxiter; 39 39 -
issm/trunk/src/c/solutions/control_core.cpp
r9393 r9563 20 20 int nsteps; 21 21 double eps_cm; 22 double tolx;23 22 bool cm_gradient; 24 23 int dim; … … 54 53 femmodel->parameters->FindParam(&cm_jump,NULL,CmJumpEnum); 55 54 femmodel->parameters->FindParam(&eps_cm,EpsCmEnum); 56 femmodel->parameters->FindParam(&tolx,TolxEnum);57 55 femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum); 58 56 femmodel->parameters->FindParam(&dim,DimEnum); … … 75 73 /*Initialize some of the BrentSearch arguments: */ 76 74 optargs.femmodel=femmodel; 77 optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx ;75 optpars.xmin=0; optpars.xmax=1; 78 76 79 77 /*Start looping: */ -
issm/trunk/src/c/solutions/solutions.h
r9218 r9563 40 40 41 41 //optimization 42 int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);43 int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);44 42 int GradJSearch(double* search_vector,FemModel* femmodel,int step); 45 43 -
issm/trunk/src/m/classes/model/model.m
r9561 r9563 184 184 maxiter = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3); 185 185 cm_responses = modelfield('default',NaN,'marshall',true,'preprocess','marshallcmresponses','format','DoubleMat','mattype',3); 186 tolx = modelfield('default',0,'marshall',true,'format','Double');187 186 optscal = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3); 188 187 eps_cm = modelfield('default',0,'marshall',true,'format','Double'); … … 191 190 cm_jump = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3); 192 191 cm_gradient = modelfield('default',0,'marshall',true,'format','Boolean'); 193 epsvel = modelfield('default',0,'marshall',true,'format','Double');194 meanvel = modelfield('default',0,'marshall',true,'format','Double');195 192 num_control_type = modelfield('default',0,'marshall',true,'format','Integer'); 196 193 num_cm_responses = modelfield('default',0,'marshall',true,'format','Integer'); … … 263 260 qmu_mass_flux_profiles = modelfield('default',NaN,'marshall',false); 264 261 qmu_mass_flux_segments = modelfield('default',{},'marshall',true,'format','MatArray'); 265 qmu_relax = modelfield('default',0,'marshall',false);266 262 267 263 %flaim … … 684 680 md.maxiter=20*ones(md.nsteps,1); 685 681 686 %tolerance used by the optimization algorithm687 md.tolx=10^-4;688 689 682 %the inversed parameter is updated as follows: 690 683 %new_par=old_par + optscal(n)*C*gradient with C in [0 1]; … … 709 702 %NaN if not applied 710 703 md.eps_cm=NaN; %not activated 711 712 %minimum velocity to avoid the misfit to be singular713 md.epsvel=eps;714 715 %averaged velocity used to scale the logarithmic Misfit (000 m/an)716 md.meanvel=1000/(365*24*3600);717 704 718 705 %grounding line migration: -
issm/trunk/src/m/model/display/displaycontrol.m
r9542 r9563 18 18 fielddisplay(md,'cm_responses','indicate the type of response for each optimization steps'); 19 19 fielddisplay(md,'maxiter','maximum iterations during each optimization step'); 20 fielddisplay(md,'tolx','minimum tolerance which will stop one optimization search');21 20 fielddisplay(md,'cm_jump','decrease threshold for misfit, default is 30%'); 22 21 fielddisplay(md,'cm_min','absolute minimum acceptable value of the inversed parameter'); 23 22 fielddisplay(md,'cm_max','absolute maximum acceptable value of the inversed parameter'); 24 23 fielddisplay(md,'cm_gradient','stop control method solution at gradient'); 25 fielddisplay(md,'meanvel','velocity scaling factor when evaluating relative or logarithmic misfit');26 fielddisplay(md,'epsvel','for relative fit, avoids misfit becoming infinity, for logarithmic fit, threshold for velocity');27 24 fielddisplay(md,'plot','visualization of the results of each iteration yes -> 1 no -> 0. Default is 1'); 28 25 disp('Available responses:'); -
issm/trunk/src/m/model/ismodelselfconsistent.m
r9542 r9563 84 84 fields={'numberofelements','numberofnodes','x','y','z','drag_coefficient','drag_p','drag_q',... 85 85 'rho_ice','rho_water','rheology_B','elementoniceshelf','surface','thickness','bed','g','lowmem','nsteps','maxiter',... 86 ' tolx','eps_res','max_nonlinear_iterations','rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface','elementconnectivity'};86 'eps_res','max_nonlinear_iterations','rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface','elementconnectivity'}; 87 87 checknan(md,fields); 88 88 %}}}} 89 89 %FIELDS >= 0 {{{1 90 90 fields={'numberofelements','numberofnodes','elements','drag_coefficient','drag_p','drag_q',... 91 'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter', 'tolx',...91 'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter',... 92 92 'lowmem','rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface'}; 93 93 checkgreater(md,fields,0); … … 95 95 %FIELDS > 0 {{{1 96 96 fields={'numberofelements','numberofnodes','elements','drag_p',... 97 'rho_ice','rho_water','rheology_B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter' ,'tolx'};97 'rho_ice','rho_water','rheology_B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter'}; 98 98 checkgreaterstrict(md,fields,0); 99 99 %}}} … … 275 275 message(['model not consistent: for qmu analysis, partitioning vector cannot go over npart, number of partition areas']); 276 276 end 277 end278 if ~md.qmu_relax,279 %if md.eps_rel>1.1*10^-5,280 % message(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);281 %end282 277 end 283 278 end -
issm/trunk/src/m/solutions/control_core.m
r9394 r9563 16 16 cm_jump=femmodel.parameters.CmJump; 17 17 eps_cm=femmodel.parameters.EpsCm; 18 tolx=femmodel.parameters.Tolx;19 18 cm_gradient=femmodel.parameters.CmGradient; 20 19 dim=femmodel.parameters.Dim; -
issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp
r8910 r9563 10 10 char* function_name=NULL; 11 11 double xmin,xmax; 12 double tolerance;13 12 double* maxiter; 14 13 OptArgs optargs; … … 33 32 FetchMatlabData(&xmin,XMIN); 34 33 FetchMatlabData(&xmax,XMAX); 35 FetchMatlabData(&tolerance,mxGetField(OPTIONS,0,"TolX"));36 34 FetchMatlabData(&maxiter,NULL,NULL,mxGetField(OPTIONS,0,"MaxIter")); 37 35 … … 45 43 optpars.xmin=xmin; 46 44 optpars.xmax=xmax; 47 optpars.tolerance=tolerance;48 45 optpars.maxiter=(int)maxiter[n_value-1]; 49 46 optpars.cm_jump=cm_jump[n_value-1];
Note:
See TracChangeset
for help on using the changeset viewer.