Changeset 9681 for issm/trunk/src/m/classes/model/model.m
- Timestamp:
- 09/08/11 08:29:34 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/m/classes/model/model.m
r9679 r9681 8 8 % {{{1 9 9 %Careful here: no other class should be used as default value this is a bug of matlab 10 cluster = modelfield('default',0,'marshall',false);11 verbose = modelfield('default',0,'marshall',true,'preprocess','marshallverbose','format','Integer');12 results = modelfield('default',0,'marshall',false);13 solver = modelfield('default',0,'marshall',false);14 debug = modelfield('default',0,'marshall',false);15 constants = modelfield('default',0,'marshall',true);16 flaim = modelfield('default',0,'marshall',false);17 surfaceforcings = modelfield('default',0,'marshall',true);18 basalforcings = modelfield('default',0,'marshall',true);19 friction = modelfield('default',0,'marshall',true);20 private = modelfield('default',0,'marshall',false);21 rifts = modelfield('default',0,'marshall',true);22 hydrology = modelfield('default',0,'marshall',true);23 settings = modelfield('default',0,'marshall',true);24 radaroverlay = modelfield('default',0,'marshall',false);25 thermal = modelfield('default',0,'marshall',true);26 miscellaneous = modelfield('default',0,'marshall',true);27 timestepping = modelfield('default',0,'marshall',true);28 groundingline = modelfield('default',0,'marshall',true);29 prognostic = modelfield('default',0,'marshall',true);30 materials = modelfield('default',0,'marshall',true);31 mask = modelfield('default',0,'marshall',true);32 qmu = modelfield('default',0,'marshall',true);10 cluster = modelfield('default',0,'marshall',false); 11 verbose = modelfield('default',0,'marshall',true,'preprocess','marshallverbose','format','Integer'); 12 results = modelfield('default',0,'marshall',false); 13 solver = modelfield('default',0,'marshall',false); 14 debug = modelfield('default',0,'marshall',false); 15 constants = modelfield('default',0,'marshall',true); 16 flaim = modelfield('default',0,'marshall',false); 17 surfaceforcings = modelfield('default',0,'marshall',true); 18 basalforcings = modelfield('default',0,'marshall',true); 19 friction = modelfield('default',0,'marshall',true); 20 private = modelfield('default',0,'marshall',false); 21 rifts = modelfield('default',0,'marshall',true); 22 hydrology = modelfield('default',0,'marshall',true); 23 settings = modelfield('default',0,'marshall',true); 24 radaroverlay = modelfield('default',0,'marshall',false); 25 thermal = modelfield('default',0,'marshall',true); 26 miscellaneous = modelfield('default',0,'marshall',true); 27 timestepping = modelfield('default',0,'marshall',true); 28 groundingline = modelfield('default',0,'marshall',true); 29 prognostic = modelfield('default',0,'marshall',true); 30 materials = modelfield('default',0,'marshall',true); 31 mask = modelfield('default',0,'marshall',true); 32 qmu = modelfield('default',0,'marshall',true); 33 33 balancethickness = modelfield('default',0,'marshall',true); 34 flowequation = modelfield('default',0,'marshall',true); 35 steadystate = modelfield('default',0,'marshall',true); 36 transient = modelfield('default',0,'marshall',true); 34 flowequation = modelfield('default',0,'marshall',true); 35 steadystate = modelfield('default',0,'marshall',true); 36 inversion = modelfield('default',0,'marshall',true); 37 transient = modelfield('default',0,'marshall',true); 37 38 diagnostic = modelfield('default',0,'marshall',true); 38 39 … … 97 98 nodeonboundary = modelfield('default',NaN,'marshall',false); 98 99 99 %Observations100 vx_obs = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);101 vy_obs = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);102 vel_obs = modelfield('default',NaN,'marshall',false);103 104 thickness_obs = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);105 106 100 %Statics parameters 107 101 connectivity = modelfield('default',0,'marshall',true,'format','Integer'); 108 102 109 %Control110 control_analysis = modelfield('default',0,'marshall',true,'format','Boolean');111 control_type = modelfield('default',NaN,'marshall',true,'preprocess','marshallcontroltype','format','DoubleMat','mattype',3);112 weights = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);113 nsteps = modelfield('default',0,'marshall',true,'format','Integer');114 maxiter = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3);115 cm_responses = modelfield('default',NaN,'marshall',true,'preprocess','marshallcmresponses','format','DoubleMat','mattype',3);116 optscal = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3);117 eps_cm = modelfield('default',0,'marshall',true,'format','Double');118 cm_min = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3);119 cm_max = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3);120 cm_jump = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3);121 cm_gradient = modelfield('default',0,'marshall',true,'format','Boolean');122 num_control_type = modelfield('default',0,'marshall',true,'format','Integer');123 num_cm_responses = modelfield('default',0,'marshall',true,'format','Integer');124 103 125 104 %Results fields … … 254 233 disp(sprintf(' Parameters:')); 255 234 disp(sprintf('%s%s%s',' Boundary conditions: type ''',inputname(1),'.bc'' to display')); 256 disp(sprintf('%s%s%s',' Observations: type ''',inputname(1),'.obs'' to display'));257 235 disp(sprintf('%s%s%s',' Materials: type ''',inputname(1),'.mat'' to display')); 258 236 disp(sprintf('%s%s%s',' Parameters: type ''',inputname(1),'.par'' to display')); … … 302 280 if isfield(structmd,'basal_melting_rate_correction'), md.basalforcings.melting_rate_correction=structmd.basal_melting_rate_correction; end 303 281 if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end 282 if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end 304 283 if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end 305 284 if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end … … 357 336 if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end 358 337 if isfield(structmd,'vertices_type'), md.flowequation.vertex_equation=structmd.vertices_type; end 359 if isfield(structmd,'eps_rel'), md. flowequation.reltol=structmd.eps_rel; end360 if isfield(structmd,'max_steadystate_iterations'), md. flowequation.maxiter=structmd.max_steadystate_iterations; end338 if isfield(structmd,'eps_rel'), md.steadystate.reltol=structmd.eps_rel; end 339 if isfield(structmd,'max_steadystate_iterations'), md.steadystate.maxiter=structmd.max_steadystate_iterations; end 361 340 if isfield(structmd,'isdiagnostic'), md.transient.isdiagnostic=structmd.isdiagnostic; end 362 341 if isfield(structmd,'isprognostic'), md.transient.isprognostic=structmd.isprognostic; end 363 342 if isfield(structmd,'isthermal'), md.transient.isthermal=structmd.isthermal; end 343 if isfield(structmd,'control_analysis'), md.inversion.iscontrol=structmd.control_analysis; end 344 if isfield(structmd,'weights'), md.inversion.cost_functions_coefficients=structmd.weights; end 345 if isfield(structmd,'nsteps'), md.inversion.nsteps=structmd.nsteps; end 346 if isfield(structmd,'maxiter_per_step'), md.inversion.maxiter_per_step=structmd.maxiter_per_step; end 347 if isfield(structmd,'cm_min'), md.inversion.min_parameters=structmd.cm_min; end 348 if isfield(structmd,'cm_max'), md.inversion.max_parameters=structmd.cm_max; end 349 if isfield(structmd,'vx_obs'), md.inversion.vx_obs=structmd.vx_obs; end 350 if isfield(structmd,'vy_obs'), md.inversion.vy_obs=structmd.vy_obs; end 351 if isfield(structmd,'vel_obs'), md.inversion.vel_obs=structmd.vel_obs; end 352 if isfield(structmd,'thickness_obs'), md.inversion.thickness_obs=structmd.thickness_obs; end 364 353 365 354 %Field changes … … 412 401 if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end 413 402 end 414 if is numeric(md.control_type),415 if (structmd.control_type==143), md. control_type={'FrictionCoefficient'}; end416 if (structmd.control_type==190), md. control_type={'RheologyBbar'}; end417 if (structmd.control_type==147), md. control_type={'Dhdt'}; end418 end 419 if is member(structmd.cm_responses(end,end),[165:170 383 388 389]),420 pos=find(structmd.cm_responses==166); md. control_type(pos)=101;421 pos=find(structmd.cm_responses==167); md. control_type(pos)=102;422 pos=find(structmd.cm_responses==168); md. control_type(pos)=103;423 pos=find(structmd.cm_responses==169); md. control_type(pos)=104;424 pos=find(structmd.cm_responses==170); md. control_type(pos)=105;425 pos=find(structmd.cm_responses==165); md. control_type(pos)=201;426 pos=find(structmd.cm_responses==389); md. control_type(pos)=501;427 pos=find(structmd.cm_responses==388); md. control_type(pos)=502;428 pos=find(structmd.cm_responses==382); md. control_type(pos)=503;403 if isfield(structmd,'control_type') & isnumeric(structmd.control_type), 404 if (structmd.control_type==143), md.inversion.control_parameters={'FrictionCoefficient'}; end 405 if (structmd.control_type==190), md.inversion.control_parameters={'RheologyBbar'}; end 406 if (structmd.control_type==147), md.inversion.control_parameters={'Thickeningrate'}; end 407 end 408 if isfield(structmd,'cost_functions') & ismember(structmd.cm_responses(end,end),[165:170 383 388 389]), 409 pos=find(structmd.cm_responses==166); md.inversion.cost_functions(pos)=101; 410 pos=find(structmd.cm_responses==167); md.inversion.cost_functions(pos)=102; 411 pos=find(structmd.cm_responses==168); md.inversion.cost_functions(pos)=103; 412 pos=find(structmd.cm_responses==169); md.inversion.cost_functions(pos)=104; 413 pos=find(structmd.cm_responses==170); md.inversion.cost_functions(pos)=105; 414 pos=find(structmd.cm_responses==165); md.inversion.cost_functions(pos)=201; 415 pos=find(structmd.cm_responses==389); md.inversion.cost_functions(pos)=501; 416 pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502; 417 pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503; 429 418 end 430 419 … … 461 450 462 451 %initialize subclasses 463 md.cluster=none; 464 md.solver=solver; 465 md.solver=addoptions(md.solver,DiagnosticVertAnalysisEnum,mumpsoptions); 466 md.results=struct(); 467 md.debug=debug; 468 md.constants=constants; 469 md.flaim=flaim; 470 md.surfaceforcings=surfaceforcings; 471 md.basalforcings=basalforcings; 472 md.friction=friction; 473 md.private=private; 474 md.rifts=rifts; 475 md.hydrology=hydrology; 476 md.settings=settings; 477 md.radaroverlay=radaroverlay; 478 md.thermal=thermal; 479 md.miscellaneous=miscellaneous; 480 md.timestepping=timestepping; 481 md.groundingline=groundingline; 482 md.prognostic=prognostic; 483 md.materials=materials; 484 md.mask=mask; 485 md.qmu=qmu; 486 md.balancethickness=balancethickness; 487 md.flowequation=flowequation; 488 md.steadystate=steadystate; 489 md.transient=transient; 452 md.cluster = none; 453 md.solver = solver; 454 md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum,mumpsoptions); 455 md.results = struct(); 456 md.debug = debug; 457 md.constants = constants; 458 md.flaim = flaim; 459 md.surfaceforcings = surfaceforcings; 460 md.basalforcings = basalforcings; 461 md.friction = friction; 462 md.private = private; 463 md.rifts = rifts; 464 md.hydrology = hydrology; 465 md.settings = settings; 466 md.radaroverlay = radaroverlay; 467 md.thermal = thermal; 468 md.miscellaneous = miscellaneous; 469 md.timestepping = timestepping; 470 md.groundingline = groundingline; 471 md.prognostic = prognostic; 472 md.materials = materials; 473 md.mask = mask; 474 md.qmu = qmu; 475 md.balancethickness = balancethickness; 476 md.flowequation = flowequation; 477 md.steadystate = steadystate; 478 md.inversion = inversion; 479 md.transient = transient; 490 480 md.diagnostic=diagnostic; 491 481 … … 508 498 md.minh=1; 509 499 510 %Control511 512 %parameter to be inferred by control methods (only513 %drag and B are supported yet)514 md.control_type={'FrictionCoefficient'};515 516 %number of steps in the control methods517 md.nsteps=20;518 519 %maximum number of iteration in the optimization algorithm for520 %each step521 md.maxiter=20*ones(md.nsteps,1);522 523 %the inversed parameter is updated as follows:524 %new_par=old_par + optscal(n)*C*gradient with C in [0 1];525 %usually the optscal must be of the order of magnitude of the526 %inversed parameter (10^8 for B, 50 for drag) and can be decreased527 %after the first iterations528 md.optscal=50*ones(md.nsteps,1);529 530 %several responses can be used:531 md.cm_responses=101*ones(md.nsteps,1);532 533 %cm_jump is used to speed up control method. When534 %misfit(1)/misfit(0) < md.cm_jump, we go directly to535 %the next step536 md.cm_jump=.7*ones(md.nsteps,1); %30 per cent decrement.537 538 %stop control solution at the gradient computation and return it?539 md.cm_gradient=0;540 541 %eps_cm is a criteria to stop the control methods.542 %if J[n]-J[n-1]/J[n] < criteria, the control run stops543 %NaN if not applied544 md.eps_cm=NaN; %not activated545 546 500 %How often to save results, default is 1 so save every step 547 501 md.output_frequency=1; 548 549 %Parallelisation parameters550 551 %cluster set as none for serial552 502 553 503 %this option can be activated to load automatically the results … … 556 506 %0 to desactivate 557 507 md.waitonlock=Inf; 558 559 508 end 560 509 %}}} … … 568 517 if(strcmp(index1.subs,'par')), displayparameters(md);return; end 569 518 if(strcmp(index1.subs,'res')), displayresults(md);return; end 570 if(strcmp(index1.subs,'obs')), displayobservations(md);return; end571 if(strcmp(index1.subs,'control')), displaycontrol(md);return; end572 519 if(strcmp(index1.subs,'parallel')), displayparallel(md);return; end 573 520 end
Note:
See TracChangeset
for help on using the changeset viewer.