multiscale_buckley_leverett_model.m

Go to the documentation of this file.

function model = multiscale_buckley_leverett_model(params)
% initialization of micromodel for buckley leverett multiscale
% problem

% Bernard Haasdonk 5.5.2010

model = [];

model = nonlin_evol_model_default;

model.mu_names = {'BL_Ul','BL_tau'};
model.mu_ranges = {[0.3,1],[1,20]};
% default values:
model.BL_Ul = 0.75;
model.BL_tau = 10;
model.T = 3*10^(-3); % Zeit, nicht Matrix
model.BL_n_x = 300;

% parameters in test-routine:
model.BL_Time=3;   % Endzeit fuer HMM-Berechnung

% Regularisierungsparameter
model.BL_epsilon = 0.00001;

model.BL_M = 2;   % water/oil viscosity ratio

%%%%%%%%%%%%%%%   makroskopische Ortsschrittweite   %%%%%%%%%%%%%%%%%%%%%%%
model.BL_DX=0.02;

%%%%%%%%%%%   Anzahl der Schritte auf mikroskopischer Skala   %%%%%%%%%%%%%
%n = n_x-2; 
model.BL_n = model.BL_n_x-2; 

%%%%%%%%%%%%%%%%%%%%%%%%   Anfangsdaten   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[U,Ur,xfront,lR,rR,X,N_X]=HMM_AB_U0_BL(DX,Ul);    
% => in init_values

%%%%%%%%%%%%%%%%%%%%%%%%%   Schrittweiten   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
compute_dt = 1;
if compute_dt
  [DT,N_T,dx,dt,LAMBDA]=HMM_CFL_BL(...
      model.BL_DX,model.BL_Time,model.BL_epsilon,model.BL_M);
  n_t=round(model.T/dt)-1;   % Anzahl mikroskopischer Zeitschritte
%  disp('model parameters should be set constant sometime');
else % variables from previous call of function...
%  DT = 
%  N_T = 
%  dx = 
%  dt = 
%  LAMBDA = 
end;

% perhaps useful:
model.BL_DT = DT;
model.BL_N_T = N_T;
model.BL_dx = dx;
model.BL_dt = dt;
model.BL_LAMBDA = LAMBDA; 

% nonlin-evol nt:
% n_t unabhaegnig von tau, Ul
model.nt = n_t+1;

% further fields required by nonlin-evol:
%model.gridtype = 'none';
model.inner_product_matrix_algorithm = @(model,model_data) ...
     eye(model.BL_n_x);
model.verbose = 10;
model.debug = 0;
model.init_values_algorithm = @BL_init_values;
model.L_E_local_ptr = @BL_L_E_local;
model.plot = @(dummy) plot(dummy);

%%% Diskretisierungsmatrix A in IR^{(nx-2)x(nx-2)}
%a = 2+dx^2/(epsilon^2*tau);
%c = epsilon^2*tau/(dx^2);
%A = c*(diag(a*ones(1,n))+diag(-ones(1,n-1),1)+diag(-ones(1,n-1),-1));
% => In detailed_data

%%% Inverse der Tridiagonalmatrix A
%[T]=INVERSE_TRIDIAG(A);
%T = inv(A);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%    % 1.Schritt: stueckweise konstante Rekonstruktion von u=RU
%    [u,x,s]=HMM_REKONSTRUKTION(U,xfront,X,dx,n_x);
%        
%    % 2.Schritt: Loesung des mikroskopischen Modells mit den Anfangsdaten u=RU
%    % auf dem mikroskopischen Gitter -> gibt Loesung u zurueck
%    
%    for iter=1:n_t    % global Loesung zum Zeitpunkt T-dt
%        [u_new,w]=HMM_MIKRO_teil(T,u,epsilon,tau,dx,dt,n_x,M);
%       if ~isequal(u_new,w*dt+u)
%         error('w*dt+u not equal new u!!');
%       end;
%       u = u_new;
%    end
%    u_sicher = u;
%            
%    i1=(15:1:20);
%    i2=(150:1:170);
%    i3=(190:1:210);
%    i4=(230:1:235);
%    index=[i1,i2,i3,i4];
%         
%    % lokale Loesung zum Zeitpunkt T
%    [u,w,index_n]=HMM_MIKRO_teil(T,u_sicher,epsilon,tau,dx,dt,n_x,M,index); 
%    X = x(index_n);
%        
%    % kontrolle/ global Loesung zum Zeitpunkt T
%    [u_kontrolle,w_k]=HMM_MIKRO_teil(T,u_sicher,epsilon,tau,dx,dt,n_x,M);
%    
%    % Fehler zwischen globaler und lokaler Loesung
%    fehler = max(abs(u_kontrolle(index_n)-u));
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  PLOT     %%

%figure  % Plot zu u-Werten
%plot(x,u_sicher,'g','LineWidth',2)
%hold on
%plot(x,u_kontrolle,'LineWidth',2)
%plot(X,u,'ro')
%title('u')
%legend('globale Lsg bei T-dt','globale Lsg bei T','lokale Lsg bei T')
%xlabel('x')
%
%figure   % Plot zu w-Werten
%plot(x,w_k*dt)
%hold on
%plot(X,w*dt,'ro')
%title('w')
%xlabel('x')

function U0 = BL_init_values(model,model_data)
%keyboard;
% Abhaegnig von Ulinks
% Unabhaegnig von Tau
% erster Eintrag Ul, rest 0!!!
U0 = zeros(1,model.BL_n_x);
U0(1)= model.BL_Ul;
%[U,Ur,xfront,lR,rR,X,N_X]=HMM_AB_U0_BL(model.BL_DX,model.BL_Ul);    
%    % 1.Schritt: stueckweise konstante Rekonstruktion von u=RU
%[U0,x,s]=HMM_REKONSTRUKTION(U,xfront,X,model.BL_dx,model.BL_n_x);

function [inc, fluxes] = BL_L_E_local(model, model_data, U, index);
if isempty(index)
  [UNew,w]=HMM_MIKRO_teil(model.T,(U(:))',...
                          model.BL_epsilon,model.BL_tau,model.BL_dx,...
                          model.dt,model.BL_n_x,model.BL_M);
  inc = w * model.dt;
else
  
end;

fluxes = [];