%TAOINVERSION class definition % % Usage: % taoinversion=taoinversion(); classdef taoinversion properties (SetAccess=public) iscontrol = 0 incomplete_adjoint = 0 control_parameters = NaN nsteps = 0 cost_functions = NaN cost_functions_coefficients = NaN min_parameters = NaN max_parameters = NaN vx_obs = NaN vy_obs = NaN vz_obs = NaN vel_obs = NaN thickness_obs = NaN end methods function obj = taoinversion(varargin) % {{{ switch nargin case 0 obj=setdefaultparameters(obj); case 1 if isa(varargin{1},'inversion'), disp('converting inversion to taoinversion'); in=varargin{1}; obj.iscontrol = in.iscontrol; obj.incomplete_adjoint = in.incomplete_adjoint; obj.control_parameters = in.control_parameters; obj.nsteps = in.nsteps; obj.cost_functions = in.cost_functions(1,:); %Keep first line only obj.cost_functions_coefficients = in.cost_functions_coefficients; obj.min_parameters = in.min_parameters; obj.max_parameters = in.max_parameters; obj.vx_obs = in.vx_obs; obj.vy_obs = in.vy_obs; obj.vz_obs = in.vz_obs; obj.vel_obs = in.vel_obs; obj.thickness_obs = in.thickness_obs; end otherwise error('constructor not supported'); end end % }}} function obj = setdefaultparameters(obj) % {{{ %default is incomplete adjoint for now obj.incomplete_adjoint=1; %parameter to be inferred by control methods (only %drag and B are supported yet) obj.control_parameters={'FrictionCoefficient'}; %number of steps in the control methods obj.nsteps=20; %several responses can be used: obj.cost_functions=101*ones(obj.nsteps,1); end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ %Early return if ~obj.iscontrol, return; end num_controls=numel(md.inversion.control_parameters); num_costfunc=size(md.inversion.cost_functions,2); md = checkfield(md,'inversion.iscontrol','values',[0 1]); md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]); md = checkfield(md,'inversion.control_parameters','cell',1,'values',... {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'}); md = checkfield(md,'inversion.nsteps','numel',1,'>=',1); md = checkfield(md,'inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]); md = checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0); md = checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]); md = checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]); if solution==BalancethicknessSolutionEnum() md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1); elseif solution==BalancethicknessSoftSolutionEnum() md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1); else md = checkfield(md,'inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1); md = checkfield(md,'inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1); end end % }}} function disp(obj) % {{{ disp(sprintf(' taoinversion parameters:')); fielddisplay(obj,'iscontrol','is inversion activated?'); fielddisplay(obj,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity'); fielddisplay(obj,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'); fielddisplay(obj,'nsteps','number of optimization searches'); fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step'); fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'); fielddisplay(obj,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'); fielddisplay(obj,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'); fielddisplay(obj,'vx_obs','observed velocity x component [m/yr]'); fielddisplay(obj,'vy_obs','observed velocity y component [m/yr]'); fielddisplay(obj,'vel_obs','observed velocity magnitude [m/yr]'); fielddisplay(obj,'thickness_obs','observed thickness [m]'); disp('Available cost functions:'); disp(' 101: SurfaceAbsVelMisfit'); disp(' 102: SurfaceRelVelMisfit'); disp(' 103: SurfaceLogVelMisfit'); disp(' 104: SurfaceLogVxVyMisfit'); disp(' 105: SurfaceAverageVelMisfit'); disp(' 201: ThicknessAbsMisfit'); disp(' 501: DragCoefficientAbsGradient'); disp(' 502: RheologyBbarAbsGradient'); disp(' 503: ThicknessAbsGradient'); end % }}} function marshall(obj,md,fid) % {{{ yts=365.0*24.0*3600.0; WriteData(fid,'object',obj,'class','inversion','fieldname','iscontrol','format','Boolean'); WriteData(fid,'enum',InversionTaoEnum(),'data',true,'format','Boolean'); if ~obj.iscontrol, return; end WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean'); WriteData(fid,'object',obj,'class','inversion','fieldname','nsteps','format','Integer'); WriteData(fid,'object',obj,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3); WriteData(fid,'object',obj,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3); WriteData(fid,'object',obj,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts); WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts); WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts); WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1); %process control parameters num_control_parameters=numel(obj.control_parameters); data=zeros(1,num_control_parameters); for i=1:num_control_parameters, data(i)=StringToEnum(obj.control_parameters{i}); end WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3); WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer'); %process cost functions num_cost_functions=size(obj.cost_functions,2); data=obj.cost_functions; pos=find(data==101); data(pos)=SurfaceAbsVelMisfitEnum(); pos=find(data==102); data(pos)=SurfaceRelVelMisfitEnum(); pos=find(data==103); data(pos)=SurfaceLogVelMisfitEnum(); pos=find(data==104); data(pos)=SurfaceLogVxVyMisfitEnum(); pos=find(data==105); data(pos)=SurfaceAverageVelMisfitEnum(); pos=find(data==201); data(pos)=ThicknessAbsMisfitEnum(); pos=find(data==501); data(pos)=DragCoefficientAbsGradientEnum(); pos=find(data==502); data(pos)=RheologyBbarAbsGradientEnum(); pos=find(data==503); data(pos)=ThicknessAbsGradientEnum(); pos=find(data==504); data(pos)=ThicknessAlongGradientEnum(); pos=find(data==505); data(pos)=ThicknessAcrossGradientEnum(); pos=find(data==506); data(pos)=BalancethicknessMisfitEnum(); WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3); WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer'); end % }}} end end