Ignore:
Timestamp:
09/08/11 08:29:34 (14 years ago)
Author:
Mathieu Morlighem
Message:

Added inversion

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/m/classes/model/model.m

    r9679 r9681  
    88                 % {{{1
    99                 %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);
    3333                 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);
    3738                 diagnostic = modelfield('default',0,'marshall',true);
    3839
     
    9798                 nodeonboundary = modelfield('default',NaN,'marshall',false);
    9899
    99                  %Observations
    100                  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 
    106100                 %Statics parameters
    107101                 connectivity             = modelfield('default',0,'marshall',true,'format','Integer');
    108102
    109                  %Control
    110                  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');
    124103
    125104                 %Results fields
     
    254233                         disp(sprintf('   Parameters:'));
    255234                         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'));
    257235                         disp(sprintf('%s%s%s','      Materials: type ''',inputname(1),'.mat'' to display'));
    258236                         disp(sprintf('%s%s%s','      Parameters: type ''',inputname(1),'.par'' to display'));
     
    302280                         if isfield(structmd,'basal_melting_rate_correction'), md.basalforcings.melting_rate_correction=structmd.basal_melting_rate_correction; end
    303281                         if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end
     282                         if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end
    304283                         if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end
    305284                         if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end
     
    357336                         if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end
    358337                         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; end
    360                          if isfield(structmd,'max_steadystate_iterations'), md.flowequation.maxiter=structmd.max_steadystate_iterations; end
     338                         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
    361340                         if isfield(structmd,'isdiagnostic'), md.transient.isdiagnostic=structmd.isdiagnostic; end
    362341                         if isfield(structmd,'isprognostic'), md.transient.isprognostic=structmd.isprognostic; end
    363342                         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
    364353
    365354                         %Field changes
     
    412401                                 if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end
    413402                         end
    414                          if isnumeric(md.control_type),
    415                                  if (structmd.control_type==143), md.control_type={'FrictionCoefficient'}; end
    416                                  if (structmd.control_type==190), md.control_type={'RheologyBbar'}; end
    417                                  if (structmd.control_type==147), md.control_type={'Dhdt'}; end
    418                          end
    419                          if ismember(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;
    429418                         end
    430419
     
    461450
    462451                         %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;
    490480                         md.diagnostic=diagnostic;
    491481
     
    508498                         md.minh=1;
    509499
    510                          %Control
    511 
    512                          %parameter to be inferred by control methods (only
    513                          %drag and B are supported yet)
    514                          md.control_type={'FrictionCoefficient'};
    515 
    516                          %number of steps in the control methods
    517                          md.nsteps=20;
    518 
    519                          %maximum number of iteration in the optimization algorithm for
    520                          %each step
    521                          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 the
    526                          %inversed parameter (10^8 for B, 50 for drag) and can be decreased
    527                          %after the first iterations
    528                          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. When
    534                          %misfit(1)/misfit(0) < md.cm_jump, we go directly to
    535                          %the next step
    536                          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 stops
    543                          %NaN if not applied
    544                          md.eps_cm=NaN; %not activated
    545 
    546500                         %How often to save results, default is 1 so save every step
    547501                         md.output_frequency=1;
    548 
    549                          %Parallelisation parameters
    550 
    551                          %cluster set as none for serial
    552502
    553503                         %this option can be activated to load automatically the results
     
    556506                         %0 to desactivate
    557507                         md.waitonlock=Inf;
    558 
    559508                 end
    560509                 %}}}
     
    568517                                 if(strcmp(index1.subs,'par')), displayparameters(md);return; end
    569518                                 if(strcmp(index1.subs,'res')), displayresults(md);return; end
    570                                  if(strcmp(index1.subs,'obs')), displayobservations(md);return; end
    571                                  if(strcmp(index1.subs,'control')), displaycontrol(md);return; end
    572519                                 if(strcmp(index1.subs,'parallel')), displayparallel(md);return; end
    573520                         end
Note: See TracChangeset for help on using the changeset viewer.