porsche_model.m

function model = porsche_model(params)
%function model = porsche_model(params)
%
% model of a Porsche Turbo
% (the geometry information is described as polygonial of k points)
%
%
% Oliver Zeeb, 28.01.11


%% ===================================
% general information
% ====================================
if nargin == 0 
  params = [];
end;

if ~isfield(params,'pdeg')
  params.pdeg = 1; % polynomial degree
end;

if ~isfield(params,'qdeg')
  params.qdeg = 2; % quadrature degree
end;

model = [];
model.verbose = 1;
model.debug = 0;
model.rb_problem_type = 'lin_stat'; 
model.compute_output_functional = 1; % turn on computation of output
model.plot_sim_data=@porsche_plot;

%% ================================
% Settings for the PDE
% --> Needed for FEM
% =================================
% if any of the following flags is set, the corresponding function
% pointers must be provided below.
model.has_source = 0; % f=0
model.has_reaction = 0;
model.has_advection = 0;
model.has_diffusivity = 1;
model.has_output_functional = 1; 
model.has_dirichlet_values = 1;
model.has_neumann_values = 1;
model.has_robin_values = 0;

model.diffusivity_tensor_coefficients=@porsche_diffusivity_tensor_coefficients;
model.diffusivity_tensor_components=@porsche_diffusivity_tensor_components;
model.diffusivity_tensor = @(glob,params) ...
    eval_affine_decomp_general(model.diffusivity_tensor_components, ... 
                               model.diffusivity_tensor_coefficients, glob,params);
%model.diffusivity_gradient = @(glob,params) zeros(size(glob,1),2);
model.dirichlet_coefficients = @(glob,params) 1;
model.dirichlet_components = @(glob,params) {zeros(size(glob,1),1)};
model.dirichlet_values = @(glob, params) ...
    eval_affine_decomp_general(model.dirichlet_components, ...
                               model.dirichlet_coefficients,glob,params);
model.neumann_coefficients = @(glob,params) 1;
model.neumann_components = @porsche_neumann_components;
model.neumann_values = @(glob,params) ...
    eval_affine_decomp_general(model.neumann_components,...
                               model.neumann_coefficients, glob,params);

model.pdeg = params.pdeg; % degree of polynomial functions
model.qdeg = params.qdeg; % quadrature degree
model.dimrange = 1; % scalar solution
model.operators = @fem_operators;

%% ==================================================================
%================== Information about the car =======================
%====================================================================
model.speed = -100; % speed of the car 
model.xrange = [0, 750]; 
model.yrange = [0, 300]; 
xmin = model.xrange(1);
xmax = model.xrange(2);
ymin = model.yrange(1);
ymax = model.yrange(2);

spline_data = load('porsche_turbo_300pts.mat');
cp = spline_data.car_spline_pts;
model.cp = cp; % cp=Car_Points
clear('spline_data');

macro_cp_indx = [1,116,138,172,199,247,276]; %indices for the macro model
macro_cp=model.cp(macro_cp_indx,:);          %points of the macro model
model.nr_macro_car_points=length(macro_cp_indx);

% macro-triangulation points on the boundary
pmacro_boundary =  ... % points on the boundary of the macro triangulation
    [xmax, ymin;         ... % punkt 8
    macro_cp(1,1), ymin; ... % punkt 9; cp(index of macro_car_point, x=1 or y=2)
    macro_cp(2,1), ymin; ... % punkt 10
    xmin, ymin;          ... % punkt 11
    xmin, macro_cp(2,2); ... % punkt 12
    xmin, macro_cp(3,2); ... % punkt 13
    xmin, ymax;          ... % punkt 14
    macro_cp(3,1), ymax; ... % punkt 15
    macro_cp(4,1), ymax; ... % punkt 16
    macro_cp(5,1), ymax; ... % punkt 17
    macro_cp(6,1), ymax; ... % punkt 18
    xmax, ymax;          ... % punkt 19
    xmax, macro_cp(6,2); ... % punkt 20
    xmax, macro_cp(7,2); ... % punkt 21
    xmax, macro_cp(1,2)];    % punkt 22
 
% all points of the macro-grid: (7 car points + 15 boundary points)
pmacro = [macro_cp; pmacro_boundary];

% Triangles of the macro-triangulation
tmacro = [9,8,1;   ... % triangle 1
          10,9,1;  ... % triangle 2
          10,1,2;  ... % triangle 3
          11,10,2; ... % triangle 4
          11,2,12; ... % triangle 5
          12,2,13; ... % triangle 6
          2,3,13;  ... % triangle 7
          13,3,14; ... % triangle 8
          3,15,14; ... % triangle 9
          3,4,15;  ... % triangle 10
          4,16,15; ... % triangle 11
          4,5,16;  ... % triangle 12
          5,17,16; ... % triangle 13
          5,18,17; ... % triangle 14
          6,18,5;  ... % triangle 15
          6,19,18; ... % triangle 16
          6,20,19; ... % triangle 17
          7,20,6;  ... % triangle 18
          7,21,20; ... % triangle 19
          22,21,7; ... % triangle 20
          1,22,7;  ... % triangle 21
          8,22,1];     % triangle 22

model.pmacro = pmacro;
model.tmacro = tmacro;
model.gridtype = 'triagrid';

%set rectangulars for defininiton of the boundary conditions:
model.bnd_rect_corner1 = [xmax-0.1, xmin-0.1, xmin-0.1, xmin-0.1, xmin+0.1; ...
                          ymin-0.1, ymin+0.1, ymin-0.1, ymax-0.1, ymin+0.1];  
model.bnd_rect_corner2 = [xmax+0.1, xmin+0.1, xmax+0.1, xmax+0.1, xmax-0.1; ...
                          ymax+0.1, ymax-0.1, ymin+0.1, xmax+0.1, ymax-0.1];
model.bnd_rect_index = [-1,-2,-2,-2,-2];
% Gamma_out, Gamma_in, Gamma_bot, Gamma_up, Gamma_int
% Dirichlet (right part), Neumann (left part), 0-Neumann (lower part), 0-Neumann (upper part), 0-Neumann (inner part)

model.use_saved_pet_for_model_data = 1; % 1 if grid data shall be loaded
                                        % from porsche_turbo_p_e_t_microgrid
model.gen_model_data = @porsche_gen_model_data;
model.gen_detailed_data = @porsche_gen_detailed_data;
model.gen_reduced_data = @porsche_gen_reduced_data;
model.detailed_simulation = @porsche_detailed_simulation;

model.output_functional = @porsche_output_functional; % detailed output-functional
model.get_dofs_from_sim_data = @(sim_data) sim_data.uh.dofs;
model.get_inner_product_matrix = @(detailed_data) ...
    detailed_data.df_info.regularized_h10_inner_product_matrix;
model.get_estimators_from_sim_data= @(sim_data) sim_data.Delta;


%% ==================================================================
%================== Reduced Basis Settings ==========================
%====================================================================
model.RB_generation_mode = 'lagrangian';
model.set_rb_in_detailed_data=@lin_stat_set_rb_in_detailed_data;
model.get_rb_size= @(model,detailed_data) size(detailed_data.RB,2);

model.decomp_mode = 0;
model.rb_simulation=@porsche_rb_simulation;
model.rb_reconstruction=@lin_stat_rb_reconstruction;

model.x1 = macro_cp(1,1);
model.x2 = macro_cp(2,1);
model.x3 = macro_cp(3,1);
model.x4 = macro_cp(4,1);
model.x5 = macro_cp(5,1);
model.x6 = macro_cp(6,1);
model.x7 = macro_cp(7,1);
model.y1 = macro_cp(1,2);
model.y2 = macro_cp(2,2);
model.y3 = macro_cp(3,2);
model.y4 = macro_cp(4,2);
model.y5 = macro_cp(5,2);
model.y6 = macro_cp(6,2);
model.y7 = macro_cp(7,2);

model.mu_names = {'x3', 'y3'}; 
%model.mu_names = {'x3', 'y3', 'x4', 'y4', 'x5', 'y5'};

model.mu_ranges = my_mu_ranges(macro_cp, model.mu_names);
model.RB_mu_list = my_RB_mu_list(model.mu_names, model.mu_ranges);

model.set_mu = @set_mu_default;
model.get_mu = @get_mu_default;
model.use_scm = 0;

model.coercivity_alpha = @porsche_coercivity_alpha;
model.operators_output = @porsche_operators_output;   % reduced output calculation

%% ==================================================================
%================== Optimization Settings ===========================
%====================================================================
model.optimization.opt_mode = 'reduced'; % or 'detailed'
model.optimization.derivatives_available = 0; % 1 if the derivatives are available after a simulation, this ist not the case for porsche_model!
model.compute_derivative_info = 1;
model.optimization.nr_constraints = 1; %= length(model.mu_names)/2;
model.optimization.init_param=set_optimization_init_param(model);
model.optimization.params_to_optimize = ones(1,length(model.mu_names));
model.optimization.objective_function = @porsche_objective_function;
model.optimization.get_Jacobian = @porsche_get_Jacobian;
model.optimization.get_Hessian = @porsche_get_Hessian_J;
% model.optimization.get_Hes_Lag = @porsche_get_Hes_Lag;
% model.optimization.get_Hessian_eq_constr = @porsche_get_Hessian_eq_constr;
% model.optimization.get_eq_constr = @porsche_get_eq_constr;
% model.optimization.get_grad_eq_constr = @porsche_get_grad_eq_constr;
model.optimization.get_Hes_Lag = @porsche_circle_x3_y3_get_Hes_Lag;
model.optimization.get_Hessian_eq_constr = @porsche_circle_x3_y3_get_Hessian_eq_constr;
model.optimization.get_eq_constr = @porsche_circle_x3_y3_get_eq_constr;
model.optimization.get_grad_eq_constr = @porsche_circle_x3_y3_get_grad_eq_constr;
model.optimization.constraint_circle_radius = 10;
model.optimization.constraint_circle_midpoint = [176.1103 ; 161.1864]; % sehr gute werte um auf mu_opt = mu_ref zu kommen

% =============================================
% resetting global functions to local functions
% needed for FEM simulation
% see elliptic_discrete_model for details
model=elliptic_discrete_model(model);
% =============================================
end


%% ==================================================================
%================== auxiliary funxtions ============================
%====================================================================
function init_param=set_optimization_init_param(model)
    init_param=zeros(1,length(model.mu_names));
    for k=1:length(model.mu_names)
        init_param(k) =  model.(model.mu_names{k});
    end
end
    

function mu_ranges = my_mu_ranges(macro_cp, mu_names)
    mu_range_all.x1 = macro_cp(1,1) + [0 0]; % x1
    mu_range_all.x2 = macro_cp(2,1) + [0 0]; % x2
    mu_range_all.x3 = macro_cp(3,1) + [-20 20]; % x3
    mu_range_all.x4 = macro_cp(4,1) + [-20 20]; % x4
    mu_range_all.x5 = macro_cp(5,1) + [-30 30]; % x5
    mu_range_all.x6 = macro_cp(6,1) + [0 0]; % x6
    mu_range_all.x7 = macro_cp(7,1) + [0 0]; % x7
    mu_range_all.y1 = macro_cp(1,2) + [0 0]; % y1
    mu_range_all.y2 = macro_cp(2,2) + [0 0]; % y2
    mu_range_all.y3 = macro_cp(3,2) + [-20 20]; % y3
    mu_range_all.y4 = macro_cp(4,2) + [-20 20]; % y4
    mu_range_all.y5 = macro_cp(5,2) + [-30 20]; % y5
    mu_range_all.y6 = macro_cp(6,2) + [0 0]; % y6
    mu_range_all.y7 = macro_cp(7,2) + [0 0]; % y7
    for k=1:length(mu_names)
        mu_ranges{k} = getfield(mu_range_all, mu_names{k});
    end
end

function RB_mu_list = my_RB_mu_list(mu_names, mu_ranges)
% old: 06.09.11
%     for k=1:length(mu_names)
%         lower_points(k) = mu_ranges{k}(1);
%         upper_points(k) = mu_ranges{k}(2);
%         middle_points(k) = (mu_ranges{k}(2) + mu_ranges{k}(1))/2;
%     end
%     RB_mu_list={lower_points, upper_points, middle_points};

      for k=1:length(mu_names)
         model_mus{k}=[mu_ranges{k}(1), (mu_ranges{k}(1)+mu_ranges{k}(2))/2, mu_ranges{k}(2)];
      end
      
      %create all combinations of this vector
      combination_mus=createCombinations(model_mus);
      
      %setting RB_mu_list
      RB_mu_list=[];
      for k=1:size(combination_mus,2)
          RB_mu_list{k} = combination_mus(:,k)';
      end
      
end