Context Navigation

source: issm/trunk/src/m/classes/taoinversion.m@ 21341

Last change on this file since 21341 was 21341, checked in by Mathieu Morlighem, 8 years ago
merged trunk-jpl and trunk for revision 21337
File size: 10.2 KB

Rev	Line
[15396]	1	%TAOINVERSION class definition
[14165]	2	%
	3	% Usage:
	4	% taoinversion=taoinversion();
	5
	6	classdef taoinversion
	7	properties (SetAccess=public)
	8	iscontrol = 0
	9	incomplete_adjoint = 0
	10	control_parameters = NaN
[17989]	11	maxsteps = 0
	12	maxiter = 0
	13	fatol = 0
	14	frtol = 0
	15	gatol = 0
	16	grtol = 0
	17	gttol = 0
	18	algorithm = ''
[14165]	19	cost_functions = NaN
	20	cost_functions_coefficients = NaN
	21	min_parameters = NaN
	22	max_parameters = NaN
	23	vx_obs = NaN
	24	vy_obs = NaN
	25	vz_obs = NaN
	26	vel_obs = NaN
	27	thickness_obs = NaN
[19105]	28	surface_obs = NaN
[14165]	29	end
	30	methods
[19105]	31	function self = extrude(self,md) % {{{
	32	self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
	33	self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
	34	self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
	35	self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
	36	if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
	37	if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
	38	if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
	39	end % }}}
	40	function self = taoinversion(varargin) % {{{
[14165]	41	switch nargin
	42	case 0
[19105]	43	self=setdefaultparameters(self);
[14170]	44	case 1
[19105]	45	self=structtoobj(taoinversion(),varargin{1});
[14165]	46	otherwise
	47	error('constructor not supported');
	48	end
	49	end % }}}
[19105]	50	function self = setdefaultparameters(self) % {{{
[14165]	51
	52	%default is incomplete adjoint for now
[19105]	53	self.incomplete_adjoint=1;
[14165]	54
	55	%parameter to be inferred by control methods (only
	56	%drag and B are supported yet)
[19105]	57	self.control_parameters={'FrictionCoefficient'};
[14165]	58
[17989]	59	%number of iterations and steps
[19105]	60	self.maxsteps=20;
	61	self.maxiter =30;
[14165]	62
[17989]	63	%default tolerances
[19105]	64	self.fatol = 0;
	65	self.frtol = 0;
	66	self.gatol = 0;
	67	self.grtol = 0;
	68	self.gttol = 1e-4;
[17989]	69
	70	%minimization algorithm
[19105]	71	PETSCMAJOR = IssmConfig('_PETSC_MAJOR_');
	72	PETSCMINOR = IssmConfig('_PETSC_MINOR_');
	73	if(PETSCMAJOR>3 \| (PETSCMAJOR==3 & PETSCMINOR>=5))
	74	self.algorithm = 'blmvm';
	75	else
	76	self.algorithm = 'tao_blmvm';
	77	end
[17989]	78
[14165]	79	%several responses can be used:
[19105]	80	self.cost_functions=101;
[14165]	81
	82	end % }}}
[19105]	83	function md = checkconsistency(self,md,solution,analyses) % {{{
[14165]	84
	85	%Early return
[19105]	86	if ~self.iscontrol, return; end
[14165]	87
[16560]	88	if ~IssmConfig('_HAVE_TAO_'),
	89	md = checkmessage(md,['TAO has not been installed, ISSM needs to be reconfigured and recompiled with TAO']);
	90	end
	91
[14165]	92	num_controls=numel(md.inversion.control_parameters);
	93	num_costfunc=size(md.inversion.cost_functions,2);
	94
[17806]	95	md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]);
	96	md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
[19105]	97	md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols());
[17989]	98	md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0);
	99	md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0);
	100	md = checkfield(md,'fieldname','inversion.fatol','numel',1,'>=',0);
	101	md = checkfield(md,'fieldname','inversion.frtol','numel',1,'>=',0);
	102	md = checkfield(md,'fieldname','inversion.gatol','numel',1,'>=',0);
	103	md = checkfield(md,'fieldname','inversion.grtol','numel',1,'>=',0);
	104	md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>=',0);
[19105]	105
	106	PETSCMAJOR = IssmConfig('_PETSC_MAJOR_');
	107	PETSCMINOR = IssmConfig('_PETSC_MINOR_');
	108	if(PETSCMAJOR>3 \| (PETSCMAJOR==3 & PETSCMINOR>=5))
	109	md = checkfield(md,'fieldname','inversion.algorithm','values',{'blmvm','cg','lmvm'});
	110	else
	111	md = checkfield(md,'fieldname','inversion.algorithm','values',{'tao_blmvm','tao_cg','tao_lmvm'});
	112	end
	113
	114	md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',supportedcostfunctions());
[17806]	115	md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
	116	md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
	117	md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
[14165]	118
[21341]	119	if strcmp(solution,'BalancethicknessSolution')
[20500]	120	md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
[21341]	121	elseif strcmp(solution,'BalancethicknessSoftSolution')
[20500]	122	md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
[14165]	123	else
[20500]	124	md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
	125	md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
[14165]	126	end
	127	end % }}}
[19105]	128	function disp(self) % {{{
[14165]	129	disp(sprintf(' taoinversion parameters:'));
[19105]	130	fielddisplay(self,'iscontrol','is inversion activated?');
	131	fielddisplay(self,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity');
	132	fielddisplay(self,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
	133	fielddisplay(self,'maxsteps','maximum number of iterations (gradient computation)');
	134	fielddisplay(self,'maxiter','maximum number of Function evaluation (forward run)');
	135	fielddisplay(self,'fatol','convergence criterion: f(X)-f(X) (X: current iteration, X: "true" solution, f: cost function)');
	136	fielddisplay(self,'frtol','convergence criterion: \|f(X)-f(X)\|/\|f(X)\|');
	137	fielddisplay(self,'gatol','convergence criterion: \|\|g(X)\|\| (g: gradient of the cost function)');
	138	fielddisplay(self,'grtol','convergence criterion: \|\|g(X)\|\|/\|f(X)\|');
	139	fielddisplay(self,'gttol','convergence criterion: \|\|g(X)\|\|/\|\|g(X0)\|\| (g(X0): gradient at initial guess X0)');
	140	fielddisplay(self,'algorithm','minimization algorithm: ''tao_blmvm'', ''tao_cg'', ''tao_lmvm''');
	141	fielddisplay(self,'cost_functions','indicate the type of response for each optimization step');
	142	fielddisplay(self,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
	143	fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex');
	144	fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex');
	145	fielddisplay(self,'vx_obs','observed velocity x component [m/yr]');
	146	fielddisplay(self,'vy_obs','observed velocity y component [m/yr]');
	147	fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]');
	148	fielddisplay(self,'thickness_obs','observed thickness [m]');
	149	fielddisplay(self,'surface_obs','observed surface elevation [m]');
[14165]	150	disp('Available cost functions:');
	151	disp(' 101: SurfaceAbsVelMisfit');
	152	disp(' 102: SurfaceRelVelMisfit');
	153	disp(' 103: SurfaceLogVelMisfit');
	154	disp(' 104: SurfaceLogVxVyMisfit');
	155	disp(' 105: SurfaceAverageVelMisfit');
	156	disp(' 201: ThicknessAbsMisfit');
	157	disp(' 501: DragCoefficientAbsGradient');
	158	disp(' 502: RheologyBbarAbsGradient');
	159	disp(' 503: ThicknessAbsGradient');
	160	end % }}}
[21341]	161	function marshall(self,prefix,md,fid) % {{{
[14165]	162
[21341]	163	yts=md.constants.yts;
[16137]	164
[21341]	165	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean');
	166	WriteData(fid,prefix,'name','md.inversion.type','data',1,'format','Integer');
[19105]	167	if ~self.iscontrol, return; end
[21341]	168	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean');
	169	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer');
	170	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer');
	171	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','fatol','format','Double');
	172	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','frtol','format','Double');
	173	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gatol','format','Double');
	174	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','grtol','format','Double');
	175	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double');
	176	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','algorithm','format','String');
	177	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
	178	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
	179	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
	180	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
	181	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
	182	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
	183	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1);
	184	WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',1);
[14165]	185
	186	%process control parameters
[19105]	187	num_control_parameters=numel(self.control_parameters);
[21341]	188	WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray');
	189	WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer');
[14165]	190
	191	%process cost functions
[19105]	192	num_cost_functions=size(self.cost_functions,2);
	193	data=marshallcostfunctions(self.cost_functions);
[21341]	194	WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray');
	195	WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer');
[14165]	196	end % }}}
	197	end
	198	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats: