Index: /issm/trunk/src/m/solutions/dakota/dakota_in_data.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_in_data.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/dakota_in_data.m (revision 32) @@ -2,7 +2,7 @@ % define the data to write the dakota .in file % -% []=dakota_in_data(md,filei,package) +% []=dakota_in_data(variables,responses,dmeth,dparams,filei,package,varargin) % -function []=dakota_in_data(md,filei,package) +function []=dakota_in_data(dmeth,variables,responses,dparams,filei,package,varargin) if ~nargin @@ -11,13 +11,56 @@ end +%% parameters + +% get default set of parameters + +params=dakota_in_params(struct()); + +% merge specified parameters into default set, whether or not +% they already exist + +fnames=fieldnames(dparams); + +for i=1:numel(fnames) + if ~isfield(params,fnames{i}) + warning('dakota_in_data:unknown_param',... + 'No parameter ''%s'' in default parameter set.',... + fnames{i}); + end + params.(fnames{i})=dparams.(fnames{i}); +end + +if strcmpi(params.analysis_driver,'matlab') && ... + isempty(params.analysis_components) + [pathstr,name,ext,versn] = fileparts(filei); + params.analysis_components=fullfile(pathstr,[name '.m' versn]); +end + +% merge method parameters, though they shouldn't be in dparams + +dmeth=dmeth_params_merge(dmeth,dparams) + + %% variables -fnames=fieldnames(md.variables); +fnames=fieldnames(variables); + for i=1:length(fnames) - fhandle=str2func([class(md.variables.(fnames{i})) '.empty']); - dvar.(fnames{i})=fhandle(); - for j=1:length(md.variables.(fnames{i})) - %call setupdesign - dvar.(fnames{i})=QmuSetupDesign(md,dvar.(fnames{i}),md.variables.(fnames{i})(j)); + +% for linear constraints, just copy + + if strcmp(class(variables.(fnames{i})),'linear_inequality_constraint') || ... + strcmp(class(variables.(fnames{i})),'linear_equality_constraint' ) + dvar.(fnames{i})=variables.(fnames{i}); + +% for variables, call the setup function + + else + fhandle=str2func([class(variables.(fnames{i})) '.empty']); + dvar.(fnames{i})=fhandle(); + for j=1:length(variables.(fnames{i})) + %call setupdesign + dvar.(fnames{i})=QmuSetupDesign(dvar.(fnames{i}),variables.(fnames{i})(j),params,varargin{:}); + end end end @@ -25,56 +68,25 @@ %% responses -fnames=fieldnames(md.responses); +fnames=fieldnames(responses); + for i=1:length(fnames) - fhandle=str2func([class(md.responses.(fnames{i})) '.empty']); - dresp.(fnames{i})=fhandle(); - for j=1:length(md.responses.(fnames{i})) - dresp.(fnames{i})(j)=md.responses.(fnames{i})(j); - end +% fhandle=str2func([class(responses.(fnames{i})) '.empty']); +% dresp.(fnames{i})=fhandle(); +% for j=1:length(responses.(fnames{i})) +% dresp.(fnames{i})(j)=responses.(fnames{i})(j); +% end + +% currently all response types can just be copied + + dresp.(fnames{i})=responses.(fnames{i}); end -%% run parameters - -params=dakota_in_params(struct()); - -params.evaluation_concurrency=md.evaluation_concurrency; -params.npart=md.npart; -disp(md.method); - -if strncmpi(md.method,'nond_l',6), - params.method ='nond_local_reliability'; - params.interval_type =md.interval_type; - params.fd_gradient_step_size=md.fd_gradient_step_size; -elseif strncmpi(md.method,'nond_s',6), - params.method ='nond_sampling'; - params.seed =md.seed; - params.samples =md.samples; - params.sample_type =md.sample_type; -elseif strncmpi(md.method,'conmin_f',8), - params.method ='conmin_frcg'; -elseif strncmpi(md.method,'conmin_m',8), - params.method ='conmin_mfd'; -else - error('dakota_in_data: error message, method not supported yet!'); -end - -if strcmpi(md.analysis_driver,'matlab') - params.analysis_driver='matlab'; - if ~isempty(md.analysis_components) - params.analysis_components=md.analysis_components; - else - params.analysis_components=filei; - end -elseif ~isempty(md.analysis_driver) - params.analysis_driver=md.analysis_driver; -else - md.analysis_driver=params.analysis_driver; -end +%% write files %Write in file -dakota_in_write(params.method,dvar,dresp,filei,params); +dakota_in_write(dmeth.method,dmeth,dvar,dresp,params,filei,varargin{:}); %Write m file -dakota_m_write(md,params.method,dvar,dresp,filei,params,package); +dakota_m_write(dmeth.method,dmeth,dvar,dresp,params,filei,package,varargin{:}); end Index: /issm/trunk/src/m/solutions/dakota/dakota_in_params.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_in_params.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/dakota_in_params.m (revision 32) @@ -5,6 +5,4 @@ % function [params]=dakota_in_params(params) - -global ISSM_DIR; if ~nargin @@ -34,51 +32,7 @@ %% method section -% if isfield(params,'method') -% method=params.method; -% end +% all method parameters are in the dakota_method class -% method-independent - -if ~isfield(params,'output') - params.output=false; -end -if ~isfield(params,'max_iterations') - params.max_iterations=100; -end -if ~isfield(params,'max_function_evaluations') - params.max_function_evaluations=100; -end -if ~isfield(params,'convergence_tolerance') - params.convergence_tolerance=1.e-4; -end - -% nondeterministic methods - -if ~isfield(params,'distribution') - params.distribution='cumulative'; -end -if ~isfield(params,'compute') - params.compute='probabilities'; -end - -% nond_sampling - -if ~isfield(params,'seed') - params.seed=1234; -end -if ~isfield(params,'samples') - params.samples=100; -end -if ~isfield(params,'sample_type') - params.sample_type='lhs'; -end -if ~isfield(params,'all_variables') - params.all_variables=false; -end -if ~isfield(params,'variance_based_decomp') - params.variance_based_decomp=false; -end - -% nond_local_reliability +%% model section %% interface section @@ -91,5 +45,8 @@ end if ~isfield(params,'analysis_driver') - params.analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh']; + params.analysis_driver=''; +end +if ~isfield(params,'analysis_components') + params.analysis_components=''; end if ~isfield(params,'parameters_file') @@ -105,14 +62,7 @@ params.file_save=true; end -if ~isfield(params,'analysis_components') - params.analysis_components=''; -end %% responses section -% gradient specifications dependent on method -% if ~isfield(params,'no_gradients') -% params.no_gradients=false; -% end if ~isfield(params,'numerical_gradients') params.numerical_gradients=true; @@ -140,7 +90,11 @@ params.id_numerical_gradients=false; end -% if ~isfield(params,'no_hessians') -% params.no_hessians=false; -% end +% hessians not fully implemented +if ~isfield(params,'numerical_hessians') + params.numerical_gradients=true; +end +if ~isfield(params,'hessian_gradient_step_size') + params.fd_gradient_step_size=0.001; +end end Index: /issm/trunk/src/m/solutions/dakota/dakota_in_parse.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_in_parse.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/dakota_in_parse.m (revision 32) @@ -108,5 +108,5 @@ ncdv=tlist; dvar.cdv=[]; - display(sprintf(' Number of Dakota %s variables=%d.',... + display(sprintf(' Number of Dakota %s variables=%d.',... 'continuous_design',ncdv)); case 'cdv_initial_point' @@ -135,5 +135,5 @@ nnuv=tlist; dvar.nuv=[]; - display(sprintf(' Number of Dakota %s variables=%d.',... + display(sprintf(' Number of Dakota %s variables=%d.',... 'normal_uncertain',nnuv)); case 'nuv_means' @@ -167,5 +167,5 @@ ncsv=tlist; dvar.csv=[]; - display(sprintf(' Number of Dakota %s variables=%d.',... + display(sprintf(' Number of Dakota %s variables=%d.',... 'continuous_state',ncsv)); case 'csv_initial_state' @@ -236,7 +236,7 @@ dresp=[]; nof =0; +nlst=0; nnic=0; nnec=0; -nlst=0; nrf =0; @@ -262,34 +262,103 @@ nof =tlist; dresp.of =[]; - display(sprintf(' Number of Dakota %s=%d.',... + display(sprintf(' Number of Dakota %s=%d.',... 'objective_functions',nof)); + case 'objective_function_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.of(i).scale_type=char(tlist(i)); + end + case 'objective_function_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.of(i).scale =tlist(i); + end + case 'multi_objective_weights' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.of(i).weight =tlist(i); + end + + case 'num_least_squares_terms' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nlst=tlist; + dresp.lst=[]; + display(sprintf(' Number of Dakota %s=%d.',... + 'least_squares_terms',nlst)); + case 'least_squares_term_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.lst(i).scale_type=char(tlist(i)); + end + case 'least_squares_term_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.lst(i).scale =tlist(i); + end + case 'least_squares_weights' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.lst(i).weight =tlist(i); + end + case 'num_nonlinear_inequality_constraints' [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); nnic=tlist; dresp.nic=[]; - display(sprintf(' Number of Dakota %s=%d.',... - 'nonlinear_inequality_constraints',nof)); + display(sprintf(' Number of Dakota %s=%d.',... + 'nonlinear_inequality_constraints',nnic)); + case 'nonlinear_inequality_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).scale_type=char(tlist(i)); + end + case 'nonlinear_inequality_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).scale =tlist(i); + end + case 'nonlinear_inequality_lower_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).lower =tlist(i); + end + case 'nonlinear_inequality_upper_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).upper =tlist(i); + end + case 'num_nonlinear_equality_constraints' [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); nnec=tlist; - display(sprintf(' Number of Dakota %s=%d.',... - 'nonlinear_equality_constraints',nof)); - case 'num_least_squares_terms' - [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); - nlst=tlist; - dresp.lst=[]; - display(sprintf(' Number of Dakota %s=%d.',... - 'least_squares_terms',nof)); + dresp.nec=[]; + display(sprintf(' Number of Dakota %s=%d.',... + 'nonlinear_equality_constraints',nnec)); + case 'nonlinear_equality_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nec(i).scale_type=char(tlist(i)); + end + case 'nonlinear_equality_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nec(i).scale =tlist(i); + end + case 'nonlinear_equality_targets' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nec(i).target =tlist(i); + end + case 'num_response_functions' [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); nrf =tlist; dresp.rf =[]; - display(sprintf(' Number of Dakota %s=%d.',... + display(sprintf(' Number of Dakota %s=%d.',... 'response_functions',nrf)); + case 'response_descriptors' [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); - for i=1:length(tlist) - dresp.rf(i).descriptor=char(tlist(i)); - end + desc=tlist; otherwise warning('responses_parse:unrec_key',... @@ -308,29 +377,64 @@ strncmpi(tokens{1}{itoken},'responses',9) -% supply default descriptors if necessary - - if isfield(dresp,'of' ) && ~isfield(dresp.of,'descriptor') - for i=1:nof - dresp.of(i).descriptor=sprintf('obj_fn_%d',i); - end - end - if isfield(dresp,'nic') && ~isfield(dresp.nic,'descriptor') - for i=1:nnic - dresp.nic(i).descriptor=sprintf('nln_ineq_con_%d',i); - end - end - if isfield(dresp,'nec') && ~isfield(dresp.nec,'descriptor') - for i=1:nnec - dresp.nec(i).descriptor=sprintf('nln_eq_con_%d',i); - end - end - if isfield(dresp,'lst') && ~isfield(dresp.lst,'descriptor') - for i=1:nlst - dresp.lst(i).descriptor=sprintf('least_sq_term_%d',i); - end - end - if isfield(dresp,'rf' ) && ~isfield(dresp.rf,'descriptor') - for i=1:nrf - dresp.rf(i).descriptor=sprintf('response_fn_%d',i); +% assign specified or supply default descriptors + + if exist('desc','var') + idesc=0; + if isfield(dresp,'of' ) + for i=1:nof + idesc=idesc+1; + dresp.of(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'lst') + for i=1:nlst + idesc=idesc+1; + dresp.lst(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'nic') + for i=1:nnic + idesc=idesc+1; + dresp.nic(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'nec') + for i=1:nnec + idesc=idesc+1; + dresp.nec(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'rf' ) + for i=1:nrf + idesc=idesc+1; + dresp.rf(i).descriptor=char(desc(idesc)); + end + end + + else + if isfield(dresp,'of' ) + for i=1:nof + dresp.of(i).descriptor=sprintf('obj_fn_%d',i); + end + end + if isfield(dresp,'lst') + for i=1:nlst + dresp.lst(i).descriptor=sprintf('least_sq_term_%d',i); + end + end + if isfield(dresp,'nic') + for i=1:nnic + dresp.nic(i).descriptor=sprintf('nln_ineq_con_%d',i); + end + end + if isfield(dresp,'nec') + for i=1:nnec + dresp.nec(i).descriptor=sprintf('nln_eq_con_%d',i); + end + end + if isfield(dresp,'rf' ) + for i=1:nrf + dresp.rf(i).descriptor=sprintf('response_fn_%d',i); + end end end @@ -413,5 +517,5 @@ % assemble the lists by response function - if exist('nrespl','var') + if exist('nrespl','var') && isfield(dresp,'rf') if ~exist('nresp','var') nresp(1:length(dresp.rf))=floor(length(nrespl)/length(dresp.rf)); @@ -424,5 +528,5 @@ end - if exist('nprobl','var') + if exist('nprobl','var') && isfield(dresp,'rf') if ~exist('nprob','var') nprob(1:length(dresp.rf))=floor(length(nprobl)/length(dresp.rf)); @@ -435,5 +539,5 @@ end - if exist('nrell' ,'var') + if exist('nrell' ,'var') && isfield(dresp,'rf') if ~exist('nrel' ,'var') nrel (1:length(dresp.rf))=floor(length(nrell )/length(dresp.rf)); @@ -446,5 +550,5 @@ end - if exist('ngrell','var') + if exist('ngrell','var') && isfield(dresp,'rf') if ~exist('ngrel','var') ngrel(1:length(dresp.rf))=floor(length(ngrell)/length(dresp.rf)); Index: /issm/trunk/src/m/solutions/dakota/dakota_in_write.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_in_write.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/dakota_in_write.m (revision 32) @@ -2,7 +2,7 @@ % write a Dakota .in file. % -% []=dakota_in_write(method,dvar,dresp,filei,params) +% []=dakota_in_write(method,dmeth,dvar,dresp,params,filei,varargin) % -function []=dakota_in_write(method,dvar,dresp,filei,params) +function []=dakota_in_write(method,dmeth,dvar,dresp,params,filei,varargin) if ~nargin @@ -13,32 +13,10 @@ % process the input parameters -if ~exist('method' ,'var') || isempty(method) - method=input('Method: sampling (s) or local reliability (l)? ','s'); - if strncmpi(method,'s',1) - method='nond_sampling'; - elseif strncmpi(method,'l',1) - method='nond_local_reliability'; - end -end - -switch method -% case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'} -% case {'npsol_sqp'} - case {'conmin_frcg','conmin_mfd'} -% case {'optpp_cg','optpp_q_newton','optpp_fd_newton',... -% 'optpp_newton','optpp_pds'} -% case {'asynch_pattern_search'} -% case {'coliny_cobyla','coliny_direct','coliny_ea',... -% 'coliny_pattern_search','coliny_solis_wets'} -% case {'ncsu_direct'} -% case {'moga','soga'} -% case {'nl2sol','nlssol_sqp','optpp_g_newton'} - case {'nond_sampling'} - case {'nond_local_reliability'} -% case {'dace','fsu_quasi_mc','fsu_cvt'} -% case {'vector_parameter_study','list_parameter_study',... -% 'centered parameter_study','multidim_parameter_study'} - otherwise - error('Unrecognized method: ''%s''.',method); +if ~exist('method','var') || isempty(method) + method=input('Method? ','s'); +end + +if ~exist('dmeth' ,'var') || isempty(dmeth) + dmeth=dakota_method(method); end @@ -69,5 +47,5 @@ % write the method section -method_write(fidi,method,dresp,params); +method_write(fidi,dmeth,dresp,params); % write the model section @@ -77,5 +55,5 @@ % write the variables section -variables_write(fidi,method,dvar); +variables_write(fidi,dmeth,dvar); % write the interface section @@ -85,5 +63,5 @@ % write the responses section -responses_write(fidi,method,dresp,params); +responses_write(fidi,dmeth,dresp,params); fclose(fidi); @@ -109,42 +87,18 @@ %% function to write the method section of the file -function []=method_write(fidi,method,dresp,params) +function []=method_write(fidi,dmeth,dresp,params) display('Writing method section of Dakota input file.'); fprintf(fidi,'method,\n'); -fprintf(fidi,'\t%s\n',method); -param_write(fidi,'\t ','output',' ','\n',params); -param_write(fidi,'\t ','max_iterations',' = ','\n',params); -param_write(fidi,'\t ','max_function_evaluations',' = ','\n',params); -param_write(fidi,'\t ','convergence_tolerance',' = ','\n',params); - -switch method -% case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'} -% case {'npsol_sqp'} - case {'conmin_frcg','conmin_mfd'} -% case {'optpp_cg','optpp_q_newton','optpp_fd_newton',... -% 'optpp_newton','optpp_pds'} -% case {'asynch_pattern_search'} -% case {'coliny_cobyla','coliny_direct','coliny_ea',... -% 'coliny_pattern_search','coliny_solis_wets'} -% case {'ncsu_direct'} -% case {'moga','soga'} -% case {'nl2sol','nlssol_sqp','optpp_g_newton'} - case {'nond_sampling'} - param_write(fidi,'\t ','seed',' = ','\n',params); - param_write(fidi,'\t ','samples',' = ','\n',params); - param_write(fidi,'\t ','sample_type',' ','\n',params); - case {'nond_local_reliability'} -% case {'dace','fsu_quasi_mc','fsu_cvt'} -% case {'vector_parameter_study','list_parameter_study',... -% 'centered parameter_study','multidim_parameter_study'} -end +fprintf(fidi,'\t%s\n',dmeth.method); + +dmeth_params_write(dmeth,fidi); % write response levels -if isfield(dresp,'rf') +if strcmp(dmeth.type,'nond') && isfield(dresp,'rf') param_write(fidi,'\t ','distribution',' ','\n',params); - [respl,probl,rell,grell]=dresp_levels(dresp.rf); + [respl,probl,rell,grell]=prop_levels(dresp.rf); if ~isempty(respl) rlev_write(fidi,'response_levels',respl); @@ -198,5 +152,5 @@ %% function to write the variables section of the file -function []=variables_write(fidi,method,dvar) +function []=variables_write(fidi,dmeth,dvar) display('Writing variables section of Dakota input file.'); @@ -206,25 +160,9 @@ % variables vary by method -switch method -% case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'} -% case {'npsol_sqp'} - case {'conmin_frcg','conmin_mfd'} - vsets_write(fidi,dvar,'cdv','csv'); -% case {'optpp_cg','optpp_q_newton','optpp_fd_newton',... -% 'optpp_newton','optpp_pds'} -% case {'asynch_pattern_search'} -% case {'coliny_cobyla','coliny_direct','coliny_ea',... -% 'coliny_pattern_search','coliny_solis_wets'} -% case {'ncsu_direct'} -% case {'moga','soga'} -% case {'nl2sol','nlssol_sqp','optpp_g_newton'} - case {'nond_sampling'} - vsets_write(fidi,dvar,'nuv','csv'); - case {'nond_local_reliability'} - vsets_write(fidi,dvar,'nuv','csv'); -% case {'dace','fsu_quasi_mc','fsu_cvt'} -% case {'vector_parameter_study','list_parameter_study',... -% 'centered parameter_study','multidim_parameter_study'} -end +vsets_write(fidi,dvar,dmeth.variables); + +% linear constraints vary by method + +lcsets_write(fidi,dvar,dmeth.lcspec); fprintf(fidi,'\n'); @@ -234,21 +172,21 @@ %% function to write variable sets -function []=vsets_write(fidi,dvar,varargin) - -for i=1:length(varargin) - switch(varargin{i}) +function []=vsets_write(fidi,dvar,variables) + +for i=1:length(variables) + switch(variables{i}) case 'cdv' - vstring='continuous_design'; + cstring='continuous_design'; case 'nuv' - vstring='normal_uncertain'; + cstring='normal_uncertain'; case 'csv' - vstring='continuous_state'; + cstring='continuous_state'; otherwise warning('vsets_write:unrec_var',... - 'Unrecognized variable ''%s''.',varargin{i}); - end - - if isfield(dvar,varargin{i}) - vlist_write(fidi,vstring,varargin{i},dvar.(varargin{i})); + 'Unrecognized variable ''%s''.',variables{i}); + end + + if isfield(dvar,variables{i}) + vlist_write(fidi,cstring,variables{i},dvar.(variables{i})); end end @@ -258,15 +196,17 @@ %% function to write variable list -function []=vlist_write(fidi,vstring,vstring2,dvar) +function []=vlist_write(fidi,cstring,cstring2,dvar) % put variables into lists for writing -vinitpt=dvar_initpt(dvar); -vlower =dvar_lower (dvar); -vupper =dvar_upper (dvar); -vmean =dvar_mean (dvar); -vstddev=dvar_stddev(dvar); -vinitst=dvar_initst(dvar); -vdesc =dvar_desc (dvar); +pinitpt=prop_initpt(dvar); +plower =prop_lower (dvar); +pupper =prop_upper (dvar); +pmean =prop_mean (dvar); +pstddev=prop_stddev(dvar); +pinitst=prop_initst(dvar); +pstype =prop_stype (dvar); +pscale =prop_scale (dvar); +pdesc =prop_desc (dvar); % write variables @@ -275,32 +215,109 @@ disp(sprintf(' Writing %d %s variables.',length(dvar),class(dvar))); -fprintf(fidi,'\t%s = %d\n',vstring,length(dvar)); -if ~isempty(vinitpt) - fprintf(fidi,'\t %s_initial_point =\n',vstring2); - vector_write(fidi,sprintf('\t '),vinitpt,6,76); -end -if ~isempty(vlower) - fprintf(fidi,'\t %s_lower_bounds =\n',vstring2); - vector_write(fidi,sprintf('\t '),vlower ,6,76); -end -if ~isempty(vupper) - fprintf(fidi,'\t %s_upper_bounds =\n',vstring2); - vector_write(fidi,sprintf('\t '),vupper ,6,76); -end -if ~isempty(vmean) - fprintf(fidi,'\t %s_means =\n',vstring2); - vector_write(fidi,sprintf('\t '),vmean ,6,76); -end -if ~isempty(vstddev) - fprintf(fidi,'\t %s_std_deviations =\n',vstring2); - vector_write(fidi,sprintf('\t '),vstddev,6,76); -end -if ~isempty(vinitst) - fprintf(fidi,'\t %s_initial_state =\n',vstring2); - vector_write(fidi,sprintf('\t '),vinitst,6,76); -end -if ~isempty(vdesc) - fprintf(fidi,'\t %s_descriptors =\n',vstring2); - vector_write(fidi,sprintf('\t '),vdesc ,6,76); +fprintf(fidi,'\t%s = %d\n',cstring,length(dvar)); +if ~isempty(pinitpt) + fprintf(fidi,'\t %s_initial_point =\n',cstring2); + vector_write(fidi,sprintf('\t '),pinitpt,6,76); +end +if ~isempty(plower) + fprintf(fidi,'\t %s_lower_bounds =\n',cstring2); + vector_write(fidi,sprintf('\t '),plower ,6,76); +end +if ~isempty(pupper) + fprintf(fidi,'\t %s_upper_bounds =\n',cstring2); + vector_write(fidi,sprintf('\t '),pupper ,6,76); +end +if ~isempty(pmean) + fprintf(fidi,'\t %s_means =\n',cstring2); + vector_write(fidi,sprintf('\t '),pmean ,6,76); +end +if ~isempty(pstddev) + fprintf(fidi,'\t %s_std_deviations =\n',cstring2); + vector_write(fidi,sprintf('\t '),pstddev,6,76); +end +if ~isempty(pinitst) + fprintf(fidi,'\t %s_initial_state =\n',cstring2); + vector_write(fidi,sprintf('\t '),pinitst,6,76); +end +if ~isempty(pstype) + fprintf(fidi,'\t %s_scale_types =\n',cstring2); + vector_write(fidi,sprintf('\t '),pstype ,6,76); +end +if ~isempty(pscale) + fprintf(fidi,'\t %s_scales =\n',cstring2); + vector_write(fidi,sprintf('\t '),pscale ,6,76); +end +if ~isempty(pdesc) + fprintf(fidi,'\t %s_descriptors =\n',cstring2); + vector_write(fidi,sprintf('\t '),pdesc ,6,76); +end + +end + +%% function to write linear constraint sets + +function []=lcsets_write(fidi,dvar,lcspec) + +for i=1:length(lcspec) + switch(lcspec{i}) + case 'lic' + cstring ='linear_inequality_constraints'; + cstring2='linear_inequality'; + case 'lec' + cstring ='linear_equality_constraints'; + cstring2='linear_equality'; + otherwise + warning('lcspec_write:unrec_lcspec',... + 'Unrecognized linear constraint ''%s''.',lcspec{i}); + end + + if isfield(dvar,lcspec{i}) + lclist_write(fidi,cstring,cstring2,dvar.(lcspec{i})); + end +end + +end + +%% function to write linear constraint list + +function []=lclist_write(fidi,cstring,cstring2,dvar) + +% put linear constraints into lists for writing + +pmatrix=prop_matrix(dvar); +plower =prop_lower (dvar); +pupper =prop_upper (dvar); +ptarget=prop_target(dvar); +pstype =prop_stype (dvar); +pscale =prop_scale (dvar); + +% write linear constraints + +disp(sprintf(' Writing %d %s linear constraints.',... + length(dvar),class(dvar))); + +if ~isempty(pmatrix) + fprintf(fidi,'\t %s_matrix =\n',cstring2); + vector_write(fidi,sprintf('\t '),pmatrix,6,76); +end +if ~isempty(plower) + fprintf(fidi,'\t %s_lower_bounds =\n',cstring2); + vector_write(fidi,sprintf('\t '),plower ,6,76); +end +if ~isempty(pupper) + fprintf(fidi,'\t %s_upper_bounds =\n',cstring2); + vector_write(fidi,sprintf('\t '),pupper ,6,76); +end +if ~isempty(ptarget) + fprintf(fidi,'\t %s_targets =\n',cstring2); + vector_write(fidi,sprintf('\t '),ptarget,6,76); +end +if ~isempty(pstype) + fprintf(fidi,'\t %s_scale_types =\n',cstring2); + vector_write(fidi,sprintf('\t '),pstype ,6,76); +end +if ~isempty(pscale) + fprintf(fidi,'\t %s_scales =\n',cstring2); + vector_write(fidi,sprintf('\t '),pscale ,6,76); end @@ -341,5 +358,5 @@ %% function to write the responses section of the file -function []=responses_write(fidi,method,dresp,params) +function []=responses_write(fidi,dmeth,dresp,params) display('Writing responses section of Dakota input file.'); @@ -349,28 +366,6 @@ % functions, gradients, and hessians vary by method -switch method -% case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'} -% case {'npsol_sqp'} - case {'conmin_frcg','conmin_mfd'} - rsets_write(fidi,dresp,'of','nic','nec'); - ghspec_write(fidi,params,'grad'); -% case {'optpp_cg','optpp_q_newton','optpp_fd_newton',... -% 'optpp_newton','optpp_pds'} -% case {'asynch_pattern_search'} -% case {'coliny_cobyla','coliny_direct','coliny_ea',... -% 'coliny_pattern_search','coliny_solis_wets'} -% case {'ncsu_direct'} -% case {'moga','soga'} -% case {'nl2sol','nlssol_sqp','optpp_g_newton'} - case {'nond_sampling'} - rsets_write(fidi,dresp,'rf'); - ghspec_write(fidi,params); - case {'nond_local_reliability'} - rsets_write(fidi,dresp,'rf'); - ghspec_write(fidi,params,'grad'); -% case {'dace','fsu_quasi_mc','fsu_cvt'} -% case {'vector_parameter_study','list_parameter_study',... -% 'centered parameter_study','multidim_parameter_study'} -end +rsets_write(fidi,dresp,dmeth.responses); +ghspec_write(fidi,params,dmeth.ghspec); fprintf(fidi,'\n'); @@ -380,32 +375,32 @@ %% function to write response sets -function []=rsets_write(fidi,dresp,varargin) +function []=rsets_write(fidi,dresp,responses) rdesc={}; -for i=1:length(varargin) - switch(varargin{i}) +for i=1:length(responses) + switch(responses{i}) case 'of' - rstring ='objective_functions'; - rstring2='objective_function'; + cstring ='objective_functions'; + cstring2='objective_function'; case 'lst' - rstring ='least_squares_terms'; - rstring2='least_squares_term'; + cstring ='least_squares_terms'; + cstring2='least_squares_term'; case 'nic' - rstring ='nonlinear_inequality_constraints'; - rstring2='nonlinear_inequality'; + cstring ='nonlinear_inequality_constraints'; + cstring2='nonlinear_inequality'; case 'nec' - rstring ='nonlinear_equality_constraints'; - rstring2='nonlinear_equality'; + cstring ='nonlinear_equality_constraints'; + cstring2='nonlinear_equality'; case 'rf' - rstring ='response_functions'; - rstring2='response_function'; + cstring ='response_functions'; + cstring2='response_function'; otherwise warning('rsets_write:unrec_resp',... - 'Unrecognized response ''%s''.',varargin{i}); + 'Unrecognized response ''%s''.',responses{i}); end - if isfield(dresp,varargin{i}) - [rdesc]=rlist_write(fidi,rstring,rstring2,dresp.(varargin{i}),rdesc); + if isfield(dresp,responses{i}) + [rdesc]=rlist_write(fidi,cstring,cstring2,dresp.(responses{i}),rdesc); end end @@ -422,18 +417,18 @@ %% function to write response list -function [rdesc]=rlist_write(fidi,rstring,rstring2,dresp,rdesc) +function [rdesc]=rlist_write(fidi,cstring,cstring2,dresp,rdesc) % put responses into lists for writing -rstype =dresp_stype (dresp); -rscale =dresp_scale (dresp); -rweight=dresp_weight(dresp); -rlower =dresp_lower (dresp); -rupper =dresp_upper (dresp); -rtarget=dresp_target(dresp); +pstype =prop_stype (dresp); +pscale =prop_scale (dresp); +pweight=prop_weight(dresp); +plower =prop_lower (dresp); +pupper =prop_upper (dresp); +ptarget=prop_target(dresp); % accumulate descriptors into list for subsequent writing -rdesc(end+1:end+numel(dresp))=dresp_desc (dresp); +rdesc(end+1:end+numel(dresp))=prop_desc (dresp); % write responses @@ -441,34 +436,34 @@ disp(sprintf(' Writing %d %s responses.',length(dresp),class(dresp))); -fprintf(fidi,'\tnum_%s = %d\n',rstring,length(dresp)); -if ~isempty(rstype) - fprintf(fidi,'\t %s_scale_types =\n',rstring2); - vector_write(fidi,sprintf('\t '),rstype ,6,76); -end -if ~isempty(rscale) - fprintf(fidi,'\t %s_scales =\n',rstring2); - vector_write(fidi,sprintf('\t '),rscale ,6,76); -end -if ~isempty(rweight) - switch rstring2 +fprintf(fidi,'\tnum_%s = %d\n',cstring,length(dresp)); +if ~isempty(pstype) + fprintf(fidi,'\t %s_scale_types =\n',cstring2); + vector_write(fidi,sprintf('\t '),pstype ,6,76); +end +if ~isempty(pscale) + fprintf(fidi,'\t %s_scales =\n',cstring2); + vector_write(fidi,sprintf('\t '),pscale ,6,76); +end +if ~isempty(pweight) + switch cstring2 case 'objective_function' - wtype='multi_objective'; + fprintf(fidi,'\t %s_weights =\n','multi_objective'); + vector_write(fidi,sprintf('\t '),pweight,6,76); case 'least_squares_term' - wtype='least_squares'; - end - fprintf(fidi,'\t %s_weights =\n',wtype); - vector_write(fidi,sprintf('\t '),rweight,6,76); -end -if ~isempty(rlower) - fprintf(fidi,'\t %s_lower_bounds =\n',rstring2); - vector_write(fidi,sprintf('\t '),rlower ,6,76); -end -if ~isempty(rupper) - fprintf(fidi,'\t %s_upper_bounds =\n',rstring2); - vector_write(fidi,sprintf('\t '),rupper ,6,76); -end -if ~isempty(rtarget) - fprintf(fidi,'\t %s_targets =\n',rstring2); - vector_write(fidi,sprintf('\t '),rtarget,6,76); + fprintf(fidi,'\t %s_weights =\n','least_squares'); + vector_write(fidi,sprintf('\t '),pweight,6,76); + end +end +if ~isempty(plower) + fprintf(fidi,'\t %s_lower_bounds =\n',cstring2); + vector_write(fidi,sprintf('\t '),plower ,6,76); +end +if ~isempty(pupper) + fprintf(fidi,'\t %s_upper_bounds =\n',cstring2); + vector_write(fidi,sprintf('\t '),pupper ,6,76); +end +if ~isempty(ptarget) + fprintf(fidi,'\t %s_targets =\n',cstring2); + vector_write(fidi,sprintf('\t '),ptarget,6,76); end @@ -477,9 +472,9 @@ %% function to write gradient and hessian specifications -function []=ghspec_write(fidi,params,varargin) +function []=ghspec_write(fidi,params,ghspec) % gradients -if findstring(varargin,'grad') +if find_string(ghspec,'grad') if params.numerical_gradients param_write(fidi,'\t','numerical_gradients','','\n',params); @@ -497,23 +492,8 @@ % hessians (no implemented methods use them yet) -if findstring(varargin,'hess') +if find_string(ghspec,'hess') error('Hessians needed by method but not provided.'); else fprintf(fidi,'\tno_hessians\n'); -end - -end - -%% function to find a string in a cell array - -function [ifound]=findstring(cells,str) - -ifound=false; - -for i=1:length(cells) - if ischar(cells{i}) && strcmp(cells{i},str) - ifound=i; - return - end end @@ -543,4 +523,5 @@ end + %% function to write a vector on multiple lines @@ -560,7 +541,11 @@ lsvec=cmax; +% transpose vector from column-wise to row-wise + +vec=vec'; + % assemble each line, flushing when necessary -for i=1:length(vec) +for i=1:numel(vec) if isnumeric(vec(i)) sitem=sprintf('%g' ,vec(i)); Index: /issm/trunk/src/m/solutions/dakota/dakota_m_write.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_m_write.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/dakota_m_write.m (revision 32) @@ -3,7 +3,7 @@ % driver. % -% []=dakota_m_write(md,method,dvar,dresp,filem,params,package) +% []=dakota_m_write(method,dmeth,dvar,dresp,params,filem,package) % -function []=dakota_m_write(md,method,dvar,dresp,filem,params,package) +function []=dakota_m_write(method,dmeth,dvar,dresp,params,filem,package,varargin) if ~nargin @@ -14,9 +14,14 @@ % process the input parameters +if ~exist('method','var') || isempty(method) + method=input('Method? ','s'); +end + +if ~exist('dmeth' ,'var') || isempty(dmeth) + dmeth=dakota_method(method); +end + if ~exist('params','var') params=[]; -end -if ~isfield(params,'npart') - params.npart=10; end @@ -42,5 +47,5 @@ % write variables into the Matlab m-file -variables_write(md,fidm,method,dvar); +variables_write(fidm,dmeth,dvar,params,varargin{:}); % write solution into the Matlab m-file @@ -50,5 +55,5 @@ % write responses into the Matlab m-file -responses_write(fidm,method,params,dresp); +responses_write(fidm,dmeth,params,dresp); % write end of the Matlab m-file @@ -84,5 +89,5 @@ end fprintf(fidm,'\tloadmodel(infile);\n\n'); -fprintf(fidm,'\tmd=qmuname(md);\n\n'); +% fprintf(fidm,'\tmd=qmuname(md);\n\n'); end @@ -90,34 +95,14 @@ %% function to write design variables into the Matlab m-file -function []=variables_write(md,fidm,method,dvar) +function []=variables_write(fidm,dmeth,dvar,params,varargin) display('Writing variables for Matlab m-file.'); -fprintf(fidm,'%% Apply the design variables.\n\n'); +fprintf(fidm,'%% Apply the variables.\n\n'); ixc=0; % variables vary by method -switch method -% case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'} -% case {'npsol_sqp'} - case {'conmin_frcg','conmin_mfd'} - ixc=vsets_write(md,fidm,ixc,dvar,'cdv','csv'); -% case {'optpp_cg','optpp_q_newton','optpp_fd_newton',... -% 'optpp_newton','optpp_pds'} -% case {'asynch_pattern_search'} -% case {'coliny_cobyla','coliny_direct','coliny_ea',... -% 'coliny_pattern_search','coliny_solis_wets'} -% case {'ncsu_direct'} -% case {'moga','soga'} -% case {'nl2sol','nlssol_sqp','optpp_g_newton'} - case {'nond_sampling'} - ixc=vsets_write(md,fidm,ixc,dvar,'nuv','csv'); - case {'nond_local_reliability'} - ixc=vsets_write(md,fidm,ixc,dvar,'nuv','csv'); -% case {'dace','fsu_quasi_mc','fsu_cvt'} -% case {'vector_parameter_study','list_parameter_study',... -% 'centered parameter_study','multidim_parameter_study'} -end +ixc=vsets_write(fidm,ixc,dvar,dmeth.variables,params,varargin{:}); end @@ -125,9 +110,9 @@ %% function to write variable sets into the Matlab m-file -function [ixc]=vsets_write(md,fidm,ixc,dvar,varargin) - -for i=1:length(varargin) - if isfield(dvar,varargin{i}) - ixc=vlist_write(md,fidm,ixc,varargin{i},dvar.(varargin{i})); +function [ixc]=vsets_write(fidm,ixc,dvar,variables,params,varargin) + +for i=1:length(variables) + if isfield(dvar,variables{i}) + ixc=vlist_write(fidm,ixc,variables{i},dvar.(variables{i}),params,varargin{:}); end end @@ -137,5 +122,5 @@ %% function to write variable list into the Matlab m-file -function [ixc]=vlist_write(md,fidm,ixc,vtype,dvar) +function [ixc]=vlist_write(fidm,ixc,vtype,dvar,params,varargin) disp(sprintf(' Writing %d %s variables.',length(dvar),class(dvar))); @@ -156,5 +141,5 @@ %now, we need a string to put in the matlab file, which will update all the design variables %for this descriptor. - [string,ixc]=QmuUpdateFunctions(md,ixc,descriptor,dvar,i); + [string,ixc]=QmuUpdateFunctions(ixc,descriptor,dvar,params,i,varargin{:}); %dump this string in the matlab file. @@ -178,9 +163,9 @@ %% function to write responses into the Matlab m-file -function []=responses_write(fidm,method,params,dresp) +function []=responses_write(fidm,dmeth,params,dresp) display('Writing responses for Matlab m-file.'); -fprintf(fidm,'%% Calculate the response functions.\n\n'); +fprintf(fidm,'%% Calculate the responses.\n\n'); ifnvals=0; @@ -191,25 +176,5 @@ % responses vary by method -switch method -% case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'} -% case {'npsol_sqp'} - case {'conmin_frcg','conmin_mfd'} - ifnvals=rsets_write(fidm,ifnvals,params,dresp,'of','nic','nec'); -% case {'optpp_cg','optpp_q_newton','optpp_fd_newton',... -% 'optpp_newton','optpp_pds'} -% case {'asynch_pattern_search'} -% case {'coliny_cobyla','coliny_direct','coliny_ea',... -% 'coliny_pattern_search','coliny_solis_wets'} -% case {'ncsu_direct'} -% case {'moga','soga'} -% case {'nl2sol','nlssol_sqp','optpp_g_newton'} - case {'nond_sampling'} - ifnvals=rsets_write(fidm,ifnvals,params,dresp,'rf'); - case {'nond_local_reliability'} - ifnvals=rsets_write(fidm,ifnvals,params,dresp,'rf'); -% case {'dace','fsu_quasi_mc','fsu_cvt'} -% case {'vector_parameter_study','list_parameter_study',... -% 'centered parameter_study','multidim_parameter_study'} -end +ifnvals=rsets_write(fidm,ifnvals,params,dresp,dmeth.responses); fprintf(fidm,'\n'); @@ -222,9 +187,9 @@ %% function to write response sets into the Matlab m-file -function [ifnvals]=rsets_write(fidm,ifnvals,params,dresp,varargin) - -for i=1:length(varargin) - if isfield(dresp,varargin{i}) - ifnvals=rlist_write(fidm,ifnvals,params,varargin{i},dresp.(varargin{i})); +function [ifnvals]=rsets_write(fidm,ifnvals,params,dresp,responses) + +for i=1:length(responses) + if isfield(dresp,responses{i}) + ifnvals=rlist_write(fidm,ifnvals,params,responses{i},dresp.(responses{i})); end end Index: /issm/trunk/src/m/solutions/dakota/qmu.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/qmu.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/qmu.m (revision 32) @@ -2,4 +2,6 @@ %INPUT function md=qmu(md,package) %Deal with coupled ISSM or Cielo/ Dakota runs, to do sensitivity analyses. + +global ISSM_DIR; %first create temporary directory in which we will work @@ -7,9 +9,9 @@ qmudir='qmu'; if exist(qmudir,'dir') - overwrite=input('Overwrite existing ''qmu'' directory? Y/N [N]: ', 's'); + overwrite=input(['Overwrite existing ''' qmudir ''' directory? Y/N [N]: '], 's'); if strncmpi(overwrite,'y',1) system(['rm -rf ' qmudir]); else - error('Existing ''qmu'' directory not overwritten'); + error('Existing ''%s'' directory not overwritten.',qmudir); end end @@ -22,8 +24,21 @@ %create m and in files for dakota -dakota_in_data(md,'qmu',package); +if ~isfield(md.qmu_params,'analysis_driver') || ... + isempty(md.qmu_params.analysis_driver) + md.qmu_params.analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh']; +end + +if (numel(md.qmu_method) > 1) + imeth=input('Which method? '); +else + imeth=1; +end +dakota_in_data(md.qmu_method(imeth),md.variables,md.responses,md.qmu_params,'qmu',package,md); +cd .. +return %call dakota system('dakota -i qmu.in -o qmu.out -e qmu.err'); +% system('export MPIRUN_NPROCS=8;mpirun -np 4 dakota -i qmu.in -o qmu.out -e qmu.err'); %parse inputs and results from dakota Index: /issm/trunk/src/m/solutions/dakota/setupdesign/QmuSetupDesign.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/QmuSetupDesign.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/QmuSetupDesign.m (revision 32) @@ -1,3 +1,3 @@ -function dvar=QmuSetupDesign(md,dvar,variables) +function dvar=QmuSetupDesign(dvar,variables,params,varargin) %get descriptor @@ -7,17 +7,33 @@ if strcmpi(descriptor,'rho_ice') - dvar=setuprhoice(md,dvar,variables); + dvar=setuprhoice(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'rho_water') + + dvar=setuprhowater(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'heat_capacity') + + dvar=setupheatcapacity(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'thermal_conductivity') + + dvar=setupthermalconductivity(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'gravity') + + dvar=setupgravity(dvar,variables,params,varargin{:}); elseif strcmpi(descriptor,'thickness') - dvar=setupthickness(md,dvar,variables); + dvar=setupthickness(dvar,variables,params,varargin{:}); elseif strcmpi(descriptor,'drag') - dvar=setupdrag(md,dvar,variables); + dvar=setupdrag(dvar,variables,params,varargin{:}); elseif strcmpi(descriptor,'riftsfriction') - dvar=setupriftsfriction(md,dvar,variables); + dvar=setupriftsfriction(dvar,variables,params,varargin{:}); else Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setupdrag.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setupdrag.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setupdrag.m (revision 32) @@ -1,3 +1,10 @@ -function dvar=setupdrag(md,dvar,variables) +function dvar=setupdrag(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end %ok, dealing with semi-discrete distributed variable. Distribute according to how many Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setupgravity.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setupgravity.m (revision 32) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setupgravity.m (revision 32) @@ -0,0 +1,5 @@ +function dvar=setupgravity(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setupheatcapacity.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setupheatcapacity.m (revision 32) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setupheatcapacity.m (revision 32) @@ -0,0 +1,5 @@ +function dvar=setupheatcapacity(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setuprhoice.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setuprhoice.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setuprhoice.m (revision 32) @@ -1,3 +1,3 @@ -function dvar=setuprhoice(md,dvar,variables) +function dvar=setuprhoice(dvar,variables,params,varargin) dvar(end+1)=variables; Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setuprhowater.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setuprhowater.m (revision 32) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setuprhowater.m (revision 32) @@ -0,0 +1,5 @@ +function dvar=setuprhowater(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setupriftsfriction.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setupriftsfriction.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setupriftsfriction.m (revision 32) @@ -1,3 +1,10 @@ -function dvar=setupriftsfriction(md,dvar,variables) +function dvar=setupriftsfriction(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end %we have several rifts. Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setupthermalconductivity.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setupthermalconductivity.m (revision 32) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setupthermalconductivity.m (revision 32) @@ -0,0 +1,5 @@ +function dvar=setupthermalconductivity(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/solutions/dakota/setupdesign/setupthickness.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/setupdesign/setupthickness.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/setupdesign/setupthickness.m (revision 32) @@ -1,3 +1,10 @@ -function dvar=setupthickness(md,dvar,variables) +function dvar=setupthickness(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end %ok, dealing with semi-discrete distributed variable. Distribute according to how many Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/QmuUpdateFunctions.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/QmuUpdateFunctions.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/QmuUpdateFunctions.m (revision 32) @@ -1,29 +1,29 @@ -function [string,ixc]=QmuUpdateFunctions(md,ixc,descriptor,dvar,i) +function [string,ixc]=QmuUpdateFunctions(ixc,descriptor,dvar,params,i,varargin) if strcmpi(descriptor,'rho_ice') - [string,ixc]=updaterho_ice(md,ixc,dvar,i); + [string,ixc]=updaterho_ice(ixc,dvar,params,i,varargin{:}); elseif strcmpi(descriptor,'rho_water') - [string,ixc]=updaterho_water(md,ixc,dvar,i); + [string,ixc]=updaterho_water(ixc,dvar,params,i,varargin{:}); elseif strcmpi(descriptor,'heatcapacity') - [string,ixc]=updateheatcapacity(md,ixc,dvar,i); + [string,ixc]=updateheatcapacity(ixc,dvar,params,i,varargin{:}); elseif strcmpi(descriptor,'thermalconductivity') - [string,ixc]=updatethermalconductivity(md,ixc,dvar,i); + [string,ixc]=updatethermalconductivity(ixc,dvar,params,i,varargin{:}); elseif strcmpi(descriptor,'gravity') - [string,ixc]=updategravity(md,ixc,dvar,i); + [string,ixc]=updategravity(ixc,dvar,params,i,varargin{:}); elseif strcmpi(descriptor,'thickness') - [string,ixc]=updatethickness(md,ixc,dvar); + [string,ixc]=updatethickness(ixc,dvar,params,varargin{:}); elseif strcmpi(descriptor,'drag') - [string,ixc]=updatedrag(md,ixc,dvar); + [string,ixc]=updatedrag(ixc,dvar,params,varargin{:}); elseif strcmpi(descriptor,'riftsfriction') - [string,ixc]=updateriftsfriction(md,ixc,dvar); + [string,ixc]=updateriftsfriction(ixc,dvar,params,varargin{:}); else - warning('QmuUpdateFunctions warning message: ','Unrecognized design variable: %s',descriptor) + warning('QmuUpdateFunctions:warning','Unrecognized design variable: %s',descriptor) end Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updatedrag.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updatedrag.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updatedrag.m (revision 32) @@ -1,3 +1,10 @@ -function [string,ixc]=updatedrag(md,ixc,dvar) +function [string,ixc]=updatedrag(ixc,dvar,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end string=sprintf('\n\t[epart,npart]=MeshPartition(md,%d);\n\n',md.npart); @@ -7,5 +14,5 @@ ixc=ixc+1; idrag =sscanf(dvar(j).descriptor(5:end),'%d'); - if strcmpi(md.analysis_driver,'matlab') + if strcmpi(params.analysis_driver,'matlab') string=[string sprintf('\tdragi(%d,1)=Dakota.xC(%d);\n',idrag,ixc)]; else Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updategravity.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updategravity.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updategravity.m (revision 32) @@ -1,6 +1,6 @@ -function [string,ixc]=updategravity(md,ixc,dvar,i) +function [string,ixc]=updategravity(ixc,dvar,params,i,varargin) ixc=ixc+1; -if strcmpi(md.analysis_driver,'matlab') +if strcmpi(params.analysis_driver,'matlab') string=sprintf('\tmd.g=Dakota.xC(%d);\n',ixc); else Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updateheatcapacity.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updateheatcapacity.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updateheatcapacity.m (revision 32) @@ -1,6 +1,6 @@ -function [string,ixc]=updateheatcapacity(md,ixc,dvar,i) +function [string,ixc]=updateheatcapacity(ixc,dvar,params,i,varargin) ixc=ixc+1; -if strcmpi(md.analysis_driver,'matlab') +if strcmpi(params.analysis_driver,'matlab') string=sprintf('\tmd.heatcapacity=Dakota.xC(%d);\n',ixc); else Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_ice.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_ice.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_ice.m (revision 32) @@ -1,6 +1,6 @@ -function [string,ixc]=updaterho_ice(md,ixc,dvar,i) +function [string,ixc]=updaterho_ice(ixc,dvar,params,i,varargin) ixc=ixc+1; -if strcmpi(md.analysis_driver,'matlab') +if strcmpi(params.analysis_driver,'matlab') string=sprintf('\tmd.rho_ice=Dakota.xC(%d);\n',ixc); else Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_water.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_water.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_water.m (revision 32) @@ -1,6 +1,6 @@ -function [string,ixc]=updaterho_water(md,ixc,dvar,i) +function [string,ixc]=updaterho_water(ixc,dvar,params,i,varargin) ixc=ixc+1; -if strcmpi(md.analysis_driver,'matlab') +if strcmpi(params.analysis_driver,'matlab') string=sprintf('\tmd.rho_water=Dakota.xC(%d);\n',ixc); else Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updateriftsfriction.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updateriftsfriction.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updateriftsfriction.m (revision 32) @@ -1,3 +1,10 @@ -function [string,ixc]=updateriftsfriction(md,ixc,dvar) +function [string,ixc]=updateriftsfriction(ixc,dvar,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end string=sprintf('\n'); @@ -6,5 +13,5 @@ ixc=ixc+1; irift=sscanf(dvar(j).descriptor(14:end),'%d'); - if strcmpi(md.analysis_driver,'matlab') + if strcmpi(params.analysis_driver,'matlab') string=[string sprintf('\tmd.rifts(%d).friction=Dakota.xC(%d);\n',irift,ixc)]; else Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updatethermalconductivity.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updatethermalconductivity.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updatethermalconductivity.m (revision 32) @@ -1,6 +1,6 @@ -function [string,ixc]=updatethermalconductivity(md,ixc,dvar,i) +function [string,ixc]=updatethermalconductivity(ixc,dvar,params,i,varargin) ixc=ixc+1; -if strcmpi(md.analysis_driver,'matlab') +if strcmpi(params.analysis_driver,'matlab') string=sprintf('\tmd.thermalconductivity=Dakota.xC(%d);\n',ixc); else Index: /issm/trunk/src/m/solutions/dakota/updatefunctions/updatethickness.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/updatefunctions/updatethickness.m (revision 31) +++ /issm/trunk/src/m/solutions/dakota/updatefunctions/updatethickness.m (revision 32) @@ -1,3 +1,10 @@ -function [string,ixc]=updatethickness(md,ixc,dvar) +function [string,ixc]=updatethickness(ixc,dvar,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end string=sprintf('\n\t[epart,npart]=MeshPartition(md,%d);\n',md.npart); @@ -7,5 +14,5 @@ ixc=ixc+1; ithick=sscanf(dvar(j).descriptor(10:end),'%d'); - if strcmpi(md.analysis_driver,'matlab') + if strcmpi(params.analysis_driver,'matlab') string=[string sprintf('\tthicki(%d,1)=Dakota.xC(%d);\n',ithick,ixc)]; else Index: /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m (revision 31) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m (revision 32) @@ -312,10 +312,10 @@ if ~isempty(md.penalties) & ~isnan(md.penalties), for i=1:size(md.penalties,1), - for j=1:(md.numlayers-1), + for j=2:size(md.penalties,2), %constrain first dof constraints(count).constraint=rgb; constraints(count).constraint.grid1=md.penalties(i,1); - constraints(count).constraint.grid2=md.penalties(i,j+1); + constraints(count).constraint.grid2=md.penalties(i,j); constraints(count).constraint.dof=1; count=count+1; @@ -324,5 +324,5 @@ constraints(count).constraint=rgb; constraints(count).constraint.grid1=md.penalties(i,1); - constraints(count).constraint.grid2=md.penalties(i,j+1); + constraints(count).constraint.grid2=md.penalties(i,j); constraints(count).constraint.dof=2; count=count+1; Index: /issm/trunk/src/m/solutions/ice/ModelProcessorPrognostic.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorPrognostic.m (revision 31) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorPrognostic.m (revision 32) @@ -116,6 +116,19 @@ loads=struct('load',cell([0,1])); -%Single point constraints: none; -constraints=struct('constraint',cell(0,0)); + +%Initialize constraints structure +constraints=struct('constraint',cell(length(find(md.gridondirichlet_prog)),1)); + +pos=find(md.gridondirichlet_prog); +count=1; +for i=1:length(find(md.gridondirichlet_prog)), + + %constrain the thickness + constraints(count).constraint=spc; + constraints(count).constraint.grid=pos(i); + constraints(count).constraint.dof=1; + constraints(count).constraint.value=md.dirichletvalues_prog(pos(i)); %in m + count=count+1; +end %Last thing, return a partitioning vector, and its transpose. Index: /issm/trunk/src/m/solutions/ice/StrainRateCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/StrainRateCompute.m (revision 31) +++ /issm/trunk/src/m/solutions/ice/StrainRateCompute.m (revision 32) @@ -34,5 +34,5 @@ strainratematrix=[strainratevector(1) strainratevector(3) - strainratevector(3) strainratevector(2)]; + strainratevector(3) strainratevector(2)]; %eigen values and vectors @@ -41,7 +41,7 @@ %Plug into global vectors strainrate1(n,:)=strainratevector; - A1(n,1)=value(1,1); A2(n,1)=value(2,2); - Vx1(n,1)=directions(1,1); Vx2(n,1)=directions(1,2); - Vy1(n,1)=directions(2,1); Vy2(n,1)=directions(2,2); + A1(n,1)=value(1,1); A2(n,1)=value(2,2); + Vx1(n,1)=directions(1,1); Vx2(n,1)=directions(1,2); + Vy1(n,1)=directions(2,1); Vy2(n,1)=directions(2,2); end end @@ -54,7 +54,7 @@ strainrate.xy=strainrate1(:,3)*(365.25*24*3600); %principal axis and values - strainrate.principalvalue2=A1; + strainrate.principalvalue2=A1*(365.25*24*3600); %strain rate in 1/a instead of 1/s strainrate.principalaxis2=[Vx1 Vy1]; - strainrate.principalvalue1=A2; + strainrate.principalvalue1=A2*(365.25*24*3600); strainrate.principalaxis1=[Vx2 Vy2]; %Norm or effective value Index: /issm/trunk/src/m/solutions/ice/icediagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic.m (revision 32) +++ /issm/trunk/src/m/solutions/ice/icediagnostic.m (revision 32) @@ -0,0 +1,89 @@ +function md=icediagnostic(md); +%ICEDIAGNOSTIC - compute the velocity field of a model +% +% this routine is a Matlab drivn solution +% +% Usage: +% u_g=icediagnostic_core_linear(m,analysis_type,varargin); +% +% See also: ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR + +%define global variables +iceglobal + +%determine if run is parallel +if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end; + +%for now, only serial support is in +if cluster, + error('icediagnostic error message: parallel support not implemented yet'); +end + +%Create fem structure (input of icediagnostic3d) +fem=struct(); +%Figure out which type of elements are present +[fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type); + +if strcmpi(md.type,'2d'), + + %First, build elements,grids,loads, etc ... for horizontal model + if fem.ishutter, + fem.m_ss=CreateFemModel(md,'surface_slope_compute'); + fem.m_dhu=CreateFemModel(md,'diagnostic_hutter'); + end + if fem.ismacayealpattyn, + fem.m_dh=CreateFemModel(md,'diagnostic_horiz'); + end + + %plug inputs in fem (existing velocity) + if ~isnan(md.vx) & ~isnan(md.vy) & fem.ismacayealpattyn, + %an input velocity is present in the model, use it to bootstrap the diagnostic core nonlinear solution. + m_dh=fem.m_dh; + fem.inputs.velocity=zeros(m_dh.gridset.gsize,1); + fem.inputs.velocity(1:6:m_dh.gridset.gsize)=md.vx/md.yts; + fem.inputs.velocity(2:6:m_dh.gridset.gsize)=md.vy/md.yts; + else + fem.inputs=struct(); + end + + %compute solution + u_g=icediagnostic2d(md,fem); + + %Load results onto model + md=Loadresults(md,fem,u_g); + +else + + %First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model + fem.m_dbv=CreateFemModel(md,'diagnostic_basevert'); + fem.m_dv=CreateFemModel(md,'diagnostic_vert'); + + if fem.ismacayealpattyn, + fem.m_dh=CreateFemModel(md,'diagnostic_horiz'); + end + if fem.ishutter, + fem.m_ss=CreateFemModel(md,'surface_slope_compute'); + fem.m_dhu=CreateFemModel(md,'diagnostic_hutter'); + end + if fem.isstokes, + fem.m_bs=CreateFemModel(md,'bed_slope_compute'); + fem.m_ds=CreateFemModel(md,'diagnostic_stokes'); + end + + %plug inputs in fem (existing velocity) + if ~isnan(md.vx) & ~isnan(md.vy) & fem.ismacayealpattyn, + %an input velocity is present in the model, use it to bootstrap the diagnostic core nonlinear solution. + m_dh=fem.m_dh; + fem.inputs.velocity=zeros(m_dh.gridset.gsize,1); + fem.inputs.velocity(1:6:m_dh.gridset.gsize)=md.vx/md.yts; + fem.inputs.velocity(2:6:m_dh.gridset.gsize)=md.vy/md.yts; + else + fem.inputs=struct(); + end + + %compute solution + u_g=icediagnostic3d(md,fem); + + %Load results onto model + md=Loadresults(md,fem,u_g); +end Index: /issm/trunk/src/m/solutions/ice/icediagnostic_wrapper.m =================================================================== --- /issm/trunk/src/m/solutions/ice/icediagnostic_wrapper.m (revision 32) +++ /issm/trunk/src/m/solutions/ice/icediagnostic_wrapper.m (revision 32) @@ -0,0 +1,36 @@ +function md=icediagnostic_wrapper(md); +%ICEDIAGNOSTIC_WRAPPER - wrapper of diagnostic solution +% +% diagnostic wrapper for diagnostic matlab driven solution. +% we need a wrapper so that we can launch the diagnostic solution using Matlab parallel, +% +% Usage: +% md=icediagnostic_wrapper(md) +% +% See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR + +%start timing +t1=clock; + +%If running in parallel, launch core diagnostic routine within a parallel job. +if strcmpi(md.cluster,'yes'), + sched = findResource('scheduler', 'configuration', md.scheduler_configuration); + pjob = createParallelJob(sched); + set(pjob, 'FileDependencies', {'icediagnostic.m'}); + set(pjob, 'MaximumNumberOfWorkers', md.np); set(pjob, 'MinimumNumberOfWorkers', md.np); + t = createTask(pjob, @icediagnostic, 1, {md}); + set(t, 'CaptureCommandWindowOutput', true) + disp(' Submitting job to cluster'); + submit(pjob); + waitForState(pjob); + get(t, 'CommandWindowOutput') + get(t,'ErrorMessage') + disp(' Done running on cluster'); +else + md=icediagnostic(md); +end + +%stop timing +t2=clock; + +disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: /issm/trunk/src/m/solutions/macayeal/macayealcontrol.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/macayealcontrol.m (revision 31) +++ /issm/trunk/src/m/solutions/macayeal/macayealcontrol.m (revision 32) @@ -142,12 +142,14 @@ index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar); - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',sqrt(u.^2+v.^2)*yts,'title','Modeled velocity',... - 'data',sqrt(adjointu.^2+adjointv.^2)*yts,'title','Adjoint vectors','colorbar#all', 'on',... - 'data',sqrt(vx_obs.^2+vy_obs.^2)*yts,'title','Observed velocity',... - 'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#4',[0 100],'figure',1);pause(1); + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',sqrt(u.^2+v.^2)*yts,'title','Modeled velocity',... + 'data',sqrt(adjointu.^2+adjointv.^2)*yts,'title','Adjoint vectors','colorbar#all', 'on',... + 'data',sqrt(vx_obs.^2+vy_obs.^2)*yts,'title','Observed velocity',... + 'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#4',[0 100],'figure',1);drawnow; + end %Keep track of u and v to use in objectivefunction_C: @@ -164,10 +166,12 @@ old_direction=direction; - %visualize direction - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);pause(1); + %visualize direction + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);drawnow; + end @@ -208,9 +212,11 @@ %visualize new distribution. - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',B*B_to_p,'title',['B at iteration ' num2str(niteration)],'caxis',[200 900],'colorbar','on','figure',3);pause(1); + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',B*B_to_p,'title',['B at iteration ' num2str(niteration)],'caxis',[200 900],'colorbar','on','figure',3);drawnow; + end %evaluate new misfit. @@ -243,12 +249,14 @@ index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,drag_coeff,drag_coeff_bar); - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',sqrt(u.^2+v.^2),'title','Modeled velocity',... - 'data',sqrt(adjointu.^2+adjointv.^2),'title','Adjoint vectors','colorbar', 'all',... - 'data',sqrt(vx_obs.^2+vy_obs.^2),'title','Observed velocity',... - 'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#3',[0 100],'figure',1);pause(1); + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',sqrt(u.^2+v.^2),'title','Modeled velocity',... + 'data',sqrt(adjointu.^2+adjointv.^2),'title','Adjoint vectors','colorbar', 'all',... + 'data',sqrt(vx_obs.^2+vy_obs.^2),'title','Observed velocity',... + 'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#3',[0 100],'figure',1);drawnow; + end %Keep track of u and v to use in objectivefunction_C: @@ -268,9 +276,11 @@ %visualize direction - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);pause(1); + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);drawnow; + end @@ -319,5 +329,7 @@ %visualize new distribution. - plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);pause(1); + if md.plot==1, + plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);drawnow; + end %evaluate new misfit.