scm_infsup_demo.m

Go to the documentation of this file.

function scm_infsup_demo()
%scm_infsup_demo()
% 
% Demo of SCM (infsup constant)
% 

% David Kreplin 05.01.2015

%less parameters

%Thermal Block Model
% model = thermalblock_model_hp();
% model.gen_model_data=@lin_stat_gen_model_data;
% model_data = model.gen_model_data(model);
% params.dirichlet_gids = model_data.df_info.dirichlet_gids;


%Laminar Flow Model
params.mesh_number = 1;
model = laminar_flow_model(params);
model.has_nonlinearity = 0;
model_data = model.gen_model_data(model);
params.dirichlet_gids = model_data.bc_info.dirichlet_gids;




operator=@ (model) model.operators(model,model_data);
innerproduct_matrix = model.get_inner_product_matrix(model_data);

params.dirichlet_gids = model_data.bc_info.dirichlet_gids;
params.Ma = inf;
params.Mp = 0;
params.NumberOfParameters = 500;
params.eps=0.01;
params.max_iterations=10;
params.parallel=0;

params.type='infsup';

SCM_data = SCM(model, params, operator,innerproduct_matrix);
%alternative (see documentation): SCM_data = SCM(model, params, model_data);

rmserror=0;
maderror = 0;
innerproduct_matrix = remove_dirichlet(innerproduct_matrix,0,params.dirichlet_gids);
fprintf('Results:\ninfsup constant      upper bound       lower bound (scm infsup constant)\n');
for i = 1:20
    mu = rand_uniform(1,model.mu_ranges);
    model = set_mu(model,mu);
    model.decomp_mode=0;
    [A,~] =  model.operators(model,model_data);
    A=remove_dirichlet(A,0,params.dirichlet_gids );
    [v,d] = fast_eigs(A'*(innerproduct_matrix'\A),innerproduct_matrix,'sm');
   infsup=sqrt(min(d));
    alpha_UB =  SCMonline(mu,model,operator,SCM_data,'ub');
    alpha_LB = SCMonline(mu,model,operator,SCM_data,'lb');
    
    fprintf('%f             %f          %f\n',infsup,alpha_UB,alpha_LB);
    
    rmserror = rmserror+(infsup-alpha_LB)^2;
    maderror = maderror + abs(infsup-alpha_LB);
end
rmserror = sqrt(rmserror/20)
maderror = maderror/20



keyboard;
end



function A=remove_dirichlet(A,mode,I)
if mode == 0
    A(I,:) = [];
    A(:,I)= [];
else
    for i = 1:length(A)
        A{1,i}(I,:) = [];
        A{1,i}(:,I) = [];
    end
end
end