thermalblock_dd_model.m

Go to the documentation of this file.

function model = thermalblock_dd_model(params)
%function model = thermalblock_dd_model(params)
% Custom thermalblock model used for domain decomposition
%
% Small model for the `\mu`-dependent PDE
%
%      `- \nabla (k(x;\mu) \nabla u(x;\mu)) = h(x;\mu)`, in `\Omega`,
%                                 `u(x;\mu) = 0`,        on `\Gamma`.
%
% - `k(x;\mu)` is the block-wise constant heatcoefficient
% - `u(x;\mu)` is the solution
% - `h(x;\mu)` is the source function

% I.Maier, 19.07.2011

if nargin<1
  params = [];
end;

if ~isfield(params,'numintervals')
  params.numintervals = 16;
end;
if ~isfield(params,'pdeg')
  params.pdeg = 1;
end;
if ~isfield(params,'detailed_train_savepath')
  params.detailed_train_savepath = '/dom_dec';
end;

model = lin_stat_model_default;
model.decomp_mode = 0; % default
model.dual_mode = 0;

% parameters:
model.mu_names = {'mu_1','mu_2','mu_3','mu_f'};
model.mu_ranges = {[0.1,10],[0.1,10],[0.1,10],[0,1]};

model.RB_generation_mode = 'lagrangian';
model.detailed_train_savepath = params.detailed_train_savepath;
model.filecache_ignore_fields_in_model = {};

model.RB_mu_list = {[0.1,0.1,0.1,1],[0.1,0.1,10,1],[0.1,10,0.1,1], ...
                    [0.1,10,10,1],[10,0.1,10,1],[10,10,0.1,1],[0.1, ...
                    5.05,5.05,0],[10,0.1,0.1,1],[0.1,5.05,10,1], ...
                    [10,0.1,10,0],[10,10,0.1,0],[5.05,0.1,0.1,1], ...
                    [0.1,10,5.05,1],[0.1,10,10,0],[5.05,10,10,1], ...
                    [10,0.1,5.05,1],[0.1,10,5.05,0],[10,0.1,0.1,0], ...
                    [0.1,0.1,10,0],[0.1,10,0.1,0],[10,5.05,0.1,1], ...
                    [5.05,0.1,5.05,1],[5.05,5.05,0.1,1],[10,5.05, ...
                    0.1,0],[5.05,10,10,0],[0.1,5.05,0.1,1]};
                    
% default parameter values:
model = set_mu(model,[1,1,1,0.35]);

% data function specification:
model.has_diffusivity = 1;
model.has_advection = 0;
model.has_reaction = 0;
model.has_source = 1;
model.has_output_functional = 0;
model.has_dirichlet_values = 1;
model.has_neumann_values = 0;
model.has_robin_values = 0;
model.compute_output_functional = 1;

% source:
model.gamma_1 = 20;
model.gamma_2 = 20;
source_coefficients = @(grid,eindices,loc,params) 80*[params.mu_f; ...
                    1-params.mu_f];
source_components = @my_source_components;

model.source = @my_source_function;
model.source_function = @(grid,eindices,loc,params)...
    eval_affine_decomp_general(source_components,source_coefficients,...
    grid,eindices,loc,params);

% diffusivity_tensor:
diffusivity_tensor_coefficients = @(grid,eincices,loc,params) ...
    [params.mu_1; params.mu_2; params.mu_3; 1];
diffusivity_tensor_components = @my_diffusivity_tensor_components;

model.diffusivity_tensor = @(grid,eindices,loc,params) ...
    eval_affine_decomp_general(diffusivity_tensor_components, ...
                               diffusivity_tensor_coefficients, ...
                               grid,eindices,loc,params);

% dirichlet_values:
dirichlet_coefficients = @(grid,eindices,face_index,lloc,params) 0;
dirichlet_components = @(grid,eindices,face_index,lloc,params) ...
    { zeros(size(eindices,1),1) };

model.dirichlet_values = @(grid,eindices,face_index,lloc,params) ...
    eval_affine_decomp_general(dirichlet_components, ...
                               dirichlet_coefficients,grid, ...
                               eindices,face_index,lloc,params);

% output function
output_coefficients = @(grid,eindices,loc,params) 16;
output_components = @(grid,eindices,loc,params) ...
    { (grid.CX(eindices) < 0.25) & (grid.CY(eindices) > 0.75) };

model.output_function = @(grid,eindices,loc,params) ...
    eval_affine_decomp_general(output_components,output_coefficients, ...
    grid,eindices,loc,params);

model.operators_output = @fem_operators_output_volume_integral;

% geometry settings:
model.gridtype = 'triagrid';
model.xnumintervals = params.numintervals;
model.ynumintervals = params.numintervals;
model.xrange = [0,1];
model.yrange = [0,1];

% discretization parameters:
model.pdeg = params.pdeg;
model.qdeg = params.pdeg * 3;
model.dimrange = 1;

% for error estimation:
model.coercivity_alpha = @my_coercivity_alpha;
model.continuity_gamma = @my_continuity_gamma;

end

%%%%%%%%%%%% auxiliary functions:

function res = my_source_components(grid,eindices,loc,params)
glob = local2global(grid,eindices,loc,params);
res = { exp( -params.gamma_1 * ( (glob(:,1) - 0.5).^2 + ...
                                 (glob(:,2) - 0.5).^2 ) ), ...
        exp( -params.gamma_2 * ( (glob(:,1) - 0.875).^2 + ...
                                 (glob(:,2) - 0.875).^2 ) ) };
end

function res = my_source_function(grid,eindices,loc,params)
if params.dual_mode
    res = params.output_function(grid,eindices,loc,params);
    if ~iscell(res)
        res = -res;
    end
else
    res = params.source_function(grid,eindices,loc,params);
end
end

function res = my_diffusivity_tensor_components(grid,eindices,loc,params)
CX = grid.CX(eindices);
CY = grid.CY(eindices);
res = {[(CX<0.5) & (CY<0.5), zeros(length(CX),2), (CX<0.5) & (CY<0.5)], ...
       [(CX>0.5) & (CY<0.5), zeros(length(CX),2), (CX>0.5) & (CY<0.5)], ...
       [(CX<0.5) & (CY>0.5), zeros(length(CX),2), (CX<0.5) & (CY>0.5)], ...
       [(CX>0.5) & (CY>0.5), zeros(length(CX),2), (CX>0.5) & (CY>0.5)] ...
      };
end

function res = my_coercivity_alpha(model)
mus = model.get_mu(model);
res = min([1;mus(1:end-1)]);
end

function res = my_continuity_gamma(model)
mus = model.get_mu(model);
res = max([1;mus(1:end-1)]);
end