quadratic_poisson_model.m

Go to the documentation of this file.

function model = quadratic_poisson_model;
% quadratic poisson equation on unit interval
%
% u on [0,1] such that
%
% - laplace (k u) + m u^2 = 1
%    u(0) = u(1) = 0
%
% parameters k and m
%
% discretization with Finite differences
%
% inner product: l2-inner product, piecewise linear elements

% B. Haasdonk 6.12.2011

% for detailed_simulation
model.mu_names = {'k','m'};
model.mu_ranges ={[0.01,1],[0,10]};
model.xrange = [0,1];
model.xnumintervals = 100;
model.dx = (model.xrange(2)-model.xrange(1))/model.xnumintervals;

% for detailed simulation
model.gen_model_data  = @my_gen_model_data;
model.detailed_simulation  = @my_detailed_simulation;
model.plot_sim_data = @my_plot_sim_data;
model.set_mu = @set_mu_default;
model.b_coefficients = @my_b_coefficients;
model.b_components = @my_b_components;
model.f_coefficients = @my_f_coefficients;
model.f_components = @my_f_components;
model.a_coefficients = @my_a_coefficients;
model.au_components = @my_au_components;

% for dummy demo_rb_gui
model.gen_detailed_data  = @my_gen_detailed_data;

%%% for real rb simulation:
model.rb_reconstruction = @my_rb_reconstruction;
model.gen_reduced_data  = @my_gen_reduced_data;
model.rb_simulation  = @my_rb_simulation;

model.get_rb_size = @(model,detailed_data) model.xnumintervals-1;
model.get_mu = @get_mu_default;
model.is_stationary = 1;
% currently no reduction of N possible, dummy subset routine:
model.reduced_data_subset = @(model, reduced_data) reduced_data;
model.get_dofs_from_sim_data = @(sim_data) sim_data.u;
model.gridtype = 'onedgrid';
model.no_lines = 0;
model.axis_equal = 0;
model.inner_product_matrix = @my_inner_product_matrix;
model.RB_Mtrain_size = 1000;
model.Newton_eps = 1e-6;
model.enable_error_estimator = 0;

%%%%%%%%% auxiliary rb functions

function model_data = my_gen_model_data(model)
params.xrange = model.xrange;
params.xnumintervals = model.xnumintervals;
model_data.grid = onedgrid(params);

function sim_data = my_detailed_simulation(model, model_data)
uk = zeros(length(model_data.grid.X)-2,1);
K = model.inner_product_matrix(model_data);
b_comp = model.b_components(model);
b_coeff = model.b_coefficients(model);
b = lincomb_sequence(b_comp, b_coeff);
f_comp = model.f_components(model);
f_coeff = model.f_coefficients(model);
f = lincomb_sequence(f_comp,f_coeff);
a_coeff = model.a_coefficients(model);
% Newton loop 
errsqr = 1e10;
k = 0;
while errsqr>model.Newton_eps;
  au_comp = model.au_components(model,uk);
  au = lincomb_sequence(au_comp,a_coeff);
  DFu = 2*au + b;
  F = (uk' * (au + b))' - f;
  min_hk = DFu\F;
  uk = uk - min_hk;
  errsqr = min_hk'*K*min_hk;
  k = k+1;
end;
sim_data.u = uk;
sim_data.iterations = k;
disp(['finished detailed_simulation with ',num2str(k),' iterations'])

function p = my_plot_sim_data(model, detailed_data, sim_data, ...
                              plot_params)
uplot = [zeros(1,size(sim_data.u,2));sim_data.u;zeros(1,size(sim_data.u,2))];
p = plot(detailed_data.grid.X,uplot','Linewidth',2);
set(gca,'xlim',model.xrange);
set(gca,'ylim',[-1.5,1.5]);
if isfield(plot_params,'title')
  title(plot_params.title);
end;
if isfield(plot_params,'legend')
  legend(plot_params.legend);
end;

function detailed_data = my_gen_detailed_data(model, model_data)
% reduced basis by POD of some trajectories
detailed_data = model_data;
% random point set:
M = rand_uniform(model.RB_Mtrain_size,model.mu_ranges);
% uniform mesh:
%[M1,M2] = meshgrid();
%M = [M1(:)'; M2(:)'];
%model.mu_ranges ={[0.01,1],[0,10]};
ntrain = size(M,2);
U = zeros(length(model_data.grid.X)-2,0);
for m = 1:ntrain
  model = set_mu(model,M(:,m));
  sim_data = detailed_simulation(model,model_data);
  U = [U,sim_data.u];
end;
% orthonormalize U:
detailed_data.K = model.inner_product_matrix(model_data);
detailed_data.RB = orthonormalize_gram_schmidt(U,detailed_data.K,1e-8);
detailed_data.b_comp = model.b_components(model);
detailed_data.f_comp = model.f_components(model);
N = size(detailed_data.RB,2);
aPhiN_comp = ...
  model.au_components(model,detailed_data.RB(:,1));
Qa = length(aPhiN_comp);
detailed_data.aPhiN_comp = cell(Qa,N);
for n = 1:N
  aPhiN_comp = model.au_components(model,detailed_data.RB(:,n));
  for q = 1:Qa
    detailed_data.aPhiN_comp{q,n} = aPhiN_comp{q};
  end;
end;

function reduced_data = my_gen_reduced_data(model, detailed_data)
N = size(detailed_data.RB,2);
reduced_data = detailed_data;
reduced_data.KN = detailed_data.RB' * detailed_data.K * detailed_data.RB;
Qb = length(detailed_data.b_comp);
reduced_data.bN_comp = zeros(N,N,Qb);
for q = 1:Qb
  reduced_data.bN_comp(:,:,q) = ...
      detailed_data.RB' * ...
      detailed_data.b_comp{q} ...
      * detailed_data.RB;
end;
Qf = length(detailed_data.f_comp);
reduced_data.fN_comp = zeros(N,Qf);
for q = 1:Qf
  reduced_data.fN_comp(:,q) = detailed_data.RB' * detailed_data.f_comp{q}; 
end;
Qa = size(detailed_data.aPhiN_comp,1);
reduced_data.aN_comp = zeros(N,N,N,Qa);
for n1 = 1:N
  for q = 1:Qa
    reduced_data.aN_comp(n1,:,:,q) = ...
        detailed_data.RB' * (detailed_data.aPhiN_comp{q,n1}(:,:)) ...
        * detailed_data.RB; 
  end;
end;

function rb_sim_data = my_rb_simulation(model, reduced_data)
%rb_sim_data = detailed_simulation(model, reduced_data);
%N = model.N;
N = size(reduced_data.KN,1);
uNk = zeros(N,1);
KN = reduced_data.KN;
% newton loop 
errsqr = 1e10;
k = 0;
b_coeff = model.b_coefficients(model);
bN = reduced_data.bN_comp * b_coeff;
f_coeff = model.f_coefficients(model);
fN = reduced_data.fN_comp * f_coeff;
a_coeff = model.a_coefficients(model);
Qa = length(a_coeff);
AN = zeros(N,N,N);
for q = 1:Qa
  AN = AN + reduced_data.aN_comp(:,:,:,q) * a_coeff(q);
end;
while errsqr>model.Newton_eps;
  AuN = zeros(N,N);
  for n = 1:N
    AuN = AuN + reshape(AN(n,:,:),[N,N]) *uNk(n);
  end;
  DFuN = 2* AuN + bN;
  FN = AuN' * uNk + bN * uNk - fN;
  min_hNk = DFuN\FN;
  uNk = uNk - min_hNk;
  errsqr = min_hNk'*KN*min_hNk;
  k = k+1;
end;
rb_sim_data.uN = uNk;
rb_sim_data.iterations = k;
disp(['finished rb_simulation with ',num2str(k),' iterations'])

function rb_sim_data = my_rb_reconstruction(model, detailed_data, ...
                                            rb_sim_data)
rb_sim_data.u = detailed_data.RB * rb_sim_data.uN;

function K = my_inner_product_matrix(model_data);
% L2 inner product matrix
K = eye(length(model_data.grid.X)-2);
K = K / model_data.grid.nelements;

%%%%%%%%%%%%%%%%%%%%%% data functions:

function b_comp = my_b_components(model)
% matrix b(xi_i,xi_j) of detailed model
% single component consisting of triangular matrix
dx = model.k * (model.dx)^(-2);
d = 2* dx * ones(model.xnumintervals-1,1);
d1 = ones(model.xnumintervals-2,1)*(-dx);
b_comp = {diag(d) + diag(d1,1) + diag(d1,-1)};

function b_coeff = my_b_coefficients(model)
% matrix b(xi_i,xi_j) of detailed model
b_coeff = [model.k];

function f_comp = my_f_components(model);
% right hand side vector f(xi_i): single cell array
f_comp = {ones(model.xnumintervals-1,1)};

function f_coeff = my_f_coefficients(model)
% right hand side vector f(xi_i)
f_coeff = [1]; % always 1

function au_comp = my_au_components(model, u)
% matrix a(u,xi_i,xi_j) decomposed
au_comp = {diag(u)};

function a_coeff = my_a_coefficients(model)
% matrix a(u,xi_i,xi_j)
a_coeff = [model.m];