thermal_block.m

Go to the documentation of this file.

%step = 1; % single detailed simulation and plot
%step = 2; % affine decomposition of FEM-matrix and RHS
step = 3; % reduced problem 

switch step
  
 case 1
  params = [];
  params.B1 = 4; params.B2 = 3;
  model = thermalblock_model(params);
    
  model = model.set_mu(model,ones(params.B1*params.B2-1,1));   % 1 in each block
%  model = model.set_mu(model,rand(size(model.get_mu(model)))+0.1);   % 1 in each block
  
  % the following is parameter-independent
  model_data = gen_model_data(model);
    
  % the following is parameter-dependent
  %f = figure;
  sim_data = detailed_simulation(model,model_data);

  if model.plot_solution == 1
        f=figure;
        plot_sim_data(model,model_data,sim_data,[]);
  end
%  keyboard;
%  close(f);

%  params.B1 = 2; params.B2 = 2;
%  model = thermalblock_model(params);
%  model_data = gen_model_data(model);
%  mu = model.get_mu(model);
%  mu(:)= 1;
%  model = model.set_mu(model,mu);
%  sim_data = detailed_simulation(model,model_data);
%  plot_sim_data(model,model_data,sim_data,[]);
  
 case 2
     
    params=[];
    params.B1 = 4;
        params.B2 = 3;
    
   detailed_data = [];
   U=[];
   
   model=thermalblock_model(params);
   
   model.plot_solution = 0;
   
   % the following is parameter-independent
   model_data = gen_model_data(model);
   
   %number of simulations (NOS)
   NOS=5;
   M = rand_uniform(NOS,model.mu_ranges);
%   M(size(M,1)+1,:)=1;
   for k=1:NOS

       % the following is parameter-dependent
       %f = figure;
       
       % random parameters in model.mu_ranges!!
       model = model.set_mu(model,M(:,k));
       %model.mus = model.mus';
       sim_data = detailed_simulation(model,model_data);

       if model.plot_solution == 1
        f=figure;
        plot_sim_data(model,model_data,sim_data,[]);
       end
       
       %U(:,k)=sim_data.u_no_stiff_springs;
       U(:,k)=sim_data.u;
   end

   % Orthonormierung von U bzgl. H1 Skalarprodukt
   % .. Skalarprodukt-Matrix: G_ij = int_Omega grad(phi_i(x)) grad(phi_j(x))
   %    Hilfsroutine G = H10_scalar_product_matrix(p,t)
   % RB = PCA_fixspace(U,[],G);
   
   p=model_data.grid.p;
   e=model_data.grid.e;
   t=model_data.grid.t;
   
   c=sim_data.c;
   a=sim_data.a;
   f=sim_data.f;
   BM=sim_data.bm;
   
   Amat = sim_data.A;
   Bmat = sim_data.B;
%   Amat = sim_data.A_matlab;
%   Fmat = sim_data.F_matlab;
%   Bmat = sim_data.B_matlab;
   
   u = sim_data.u;
   ud = sim_data.ud;
   u_tmp = sim_data.u_tmp;
%   
%    K = sim_data.K_matlab;
%    M = sim_data.M_matlab;
%    F = sim_data.F_matlab;
%    
%    % Neumann:
%    Q = sim_data.Q_matlab;
%    G = sim_data.G_matlab;
%    % Dirichlet:
%    H = sim_data.H_matlab;
%    B = sim_data.B_matlab;
%    
   diri_values=find(p(2,:)==1);
   diri_values(:);
   
   ZERO=sparse(length(p), length(p));
   for i=diri_values
       ZERO(i,i)=1;
   end
   Zero_comp{1,1}=ZERO;
   
   G=H10_scalar_product_matrix(p,t,c,a,f);
   
   % check: f(x,y) = x*(1-x)*y*(1-y)
   %  H10_norm per hand und DOF-Vektor F durch Punktauswertungen
   % f = @(x,y) x.*(1-x).*y.*(1-y);
   % F = f(p(1,:),p(2,:));
   % H10_norm = F'* G * F;
   % should be identical to hand-calculation
   
   RB = PCA_fixspace(U,[],G);

% %    % => Ergebnis RB' * G * RB = eye;
% %    % check result
%    check=[];
%    check=RB'*G*RB;  %should be eye(NOS);
% %    trace_check=trace(check);    %should be NOS
% %    matrix_check=sum(sum(check)); %should be NOS
%    err  = max(max(abs(check - eye(NOS))));
%    
%    if err < 1e-6
%        disp('RB passed');
%    else 
%        disp('RB problem');
%    end

   check=[];
   check=RBcheck(RB,G,NOS);
   if check ==1
       disp('RB passed')
   else
       disp('RB Prob')
   end
       
    detailed_data.RB = RB;

   % PARAMETERUNABHAENGIG
   %       
 %   % fem-systemmatrix: komponenten zerlegung, spaerliche Matrizen
 %   % ...
 %   
 %   % glob = dreiecksmittelpunkte
 glob = get_triangle_midpoints(p,t);
 %   % {c1...c_Qc} = model.diffusivity_components_ptr(glob,model);
 [C,Cmatrix]=model.diffusivity_components_ptr(glob,model);
 %   
 %   % K1 = assema(c1,...)
 %   % ...
 %   % K_QK = assema(c_QK,...)
 % 
 %   %K_components = { K1, ... KQK};
 
 
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%  %%%adding spring constant
%  % Some scaling
%  d1=full(diag(H*H'));
%   l=find(d1==0);
%   d1(l)=Inf*ones(size(l));
%   id1=diag(sparse(1./d1));
%   % Form normal equations
%   R=H'*id1*R;
%   H=H'*id1*H;
% 
%   % Incorporate constraints (Dirichlet conditions)
%   % Find a positive spring constant of suitable size.
% 
%   np=size(p,2);
%   tscale=1/sqrt(np); % Triangle size
%   kscale=max(max(abs(K)));
%   mscale=max(max(abs(M)));
%   qscale=max(max(abs(Q)));
%   hscale=1;
% 
%   % Spring constant
%   spc=max(max(kscale*tscale^2,mscale),qscale*tscale)/(hscale*tscale^3);
% 
%  %M=F+G+spc*R;       %R



%    [K,dummy1,dummy2]=assema(p,t,sim_data.c,a,f);      
%        => K mit loescher der Dirichlet Zeilen und Spalten 
%           identisch zu Amat


 K_comp=cell(1,size(C,2));
 for i=1:size(C,2) 
     clear K K1
     [K,dummy1,dummy2]=assema(p,t,C{1,i}',a,f);
%     K_comp{1,i}=K;
%  end
    %K1=K+M+Q+spc*H;     %A 
     K1=K;
     for j=diri_values
         K1(j,:);
         K1(:,j);
     end
     K_comp{1,i}=K1;  %vorher K
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %   
 %   A_fem_components = {K_compo, M_components, Q_components};
 
 % later: M_comp and Q_comp for general elliptic problems 
 %thermalblock case M,Q are zero
 
 M_comp = {};
 Q_comp = {};
 %A_comp = cell(1,1);
 A_comp = [Zero_comp,K_comp, Q_comp, M_comp];
 
 %        
 %   % fem-rechte Seite: komponenten zerlegung
 
 %G_i = int_boundary g phi_i
 %compute with assemb
 differentneumannvalues=1;        %neumann value 1 at the bottom
 g_comp=cell(1,length(differentneumannvalues));
 for i=1:differentneumannvalues
     [dummy1,G,dummy2,dummy3]=assemb(BM,p,e);  %BM depends on differentneumannvalues
     
     %M1=F+G+spc*R;       %R
     M1=G;
     g_comp{1,i}=M1;    %M1 war vorher G
 end
 
 G_comp=[g_comp];
 %[G]=model.neumann_components_ptr(glob,model);
 
 % G wie vorher, nur in cell array packen:
 %G = ...
 
 %G_comp = {G};
 %thermalblock f is zero
 F_comp = {};
 
 B_comp = [G_comp, F_comp];
 
    % reduced_data darf keine hochdimensionalen Daten enthalten     
 detailed_data.A_components = A_comp;     
 detailed_data.B_components = B_comp;
 
 %   % check if affine decomposition works:
 
 % set mu arbitrary and multiply, different: 
 % model = model.set_mu(...)  
 % sim_data = detailed_simulation(...);

  % PARAMETERABHAENGIG
     
 theta_A = model.diffusivity_coefficients_ptr([],model);
 theta_A = [theta_A,1];
 
 theta_A = [1,theta_A];
 
 theta_G = model.neumann_value_coefficients([],model);
 theta_F = [];
 theta_B = [theta_G, theta_F];

 %   A = lincomb_sequence(A_comp, theta_A(mu));
 A = lincomb_sequence(A_comp, theta_A);
 %A = checkA(A,Amat);
 
 B = lincomb_sequence(B_comp, theta_B);
 
 u2 = A\B;   
% u2 = u2+1;

 pdeplot(p,e,t,'xydata',u,'colormap','jet','contour','on')
 title('u')
 figure,pdeplot(p,e,t,'xydata',u2,'colormap','jet','contour','on')
 title('u2')

 keyboard;
%  
%    % Goal: u2 == sim_data.u;
% if u2==sim_data.u
%      disp('Yipi');
% else
%      error('solutions not equal...');
% end;

% checkB=abs(R-Fmat);
% if find(checkB > 1e-10)              %shoud be an empty matrix
%     disp('Problem with B')   %when true, some values are above our limit
% else
%     disp('Eureka! Right hand sides are nearly equal')
% end
% 
% checkA=abs(A-Amat);
% if find(checkA > 1e-10)              %shoud be an empty matrix
%     disp('Problem with A')   %when true, some values are above our limit
% else
%     disp('Eureka! Matrices are nearly equal')
% end
% 
% checku=abs(u-u2);
% if find(checku > 1e-10)              %shoud be an empty matrix
%     disp('Problem with solution')   %when true, some values are above our limit
% else
%     disp('Eureka! Solution are nearly equal')
% end

  case 3
       
   % reduzierte Daten erzeugen: 
    for q = 1:length(A_comp)
       reduced_data.AN_components{q} = RB'*A_comp{q}*RB;
    end
   reduced_data.BN_components = {B_comp{1}'*RB};
    
   
   % reduzierte Simulation:
%    AN = linearcomb(reduced_data.AN_components, theta_A);
%    BN = linearcomb(reduced_data.BN_components, theta_B);
%    uN = AN\BN;
   
   AN = lincomb_sequence(reduced_data.AN_components, theta_A);
   BN = lincomb_sequence(reduced_data.BN_components, theta_B);
   uN = AN\BN(:);
   
   
   u3 = RB * uN;
%   % reduced output sN = ...
%   
%   if (u2==sim_data.u) && (A == sim_data.A) && (B == sim_data.B)
%     disp('Yipi');
%   else
%     error('solutions not equal...');
%   end;
%   
%   %  rb_sim_data = rb_simulation(model,reduced_data);       
%     %  in elliptic_compliant_rb_simulation
%     %  ... koeffizientenfunktionen auswerten
%     %  ... linearkombinationen bilden
%
%     
%   %model.RB_generation_mode = 'from_detailed_data';     
%   %  detailed_data = gen_detailed_data(model,model_data)    
%        
%   
%   % reduzierte basis generieren
 end;
% 
% 
% %model = thermalblock_model(params); % struktur, die Einstellungen enthaelt
% %  oder
% %model = thermalblock_model; % struktur, die Einstellungen enthaelt
% %model_data = gen_model_data(model); % Gitter anlegen und in model_data speichern
% %sim_data = detailed_simulation(model,model_data)
% 
% %pdeplot(model_data.grid.p,model_data.grid.e,model_data.grid.t,'xydata', sim_data.u);