clc; clear; close all;

Globals2D;

tic;
Problem;

N = 6; % degree of polynomial representing solution within an element (local solution representation)

StartUp2D; % computing local operators and mesh information
[mapI, vmapI, mapO, vmapO] = BuildBCMaps();

u = IC(x,y); % evaluating initial condition

alpha_flux = 1; % central flux
alpha_flux = 0; % upwind flux
preprocessing_time = toc;

tic;
% u = Timestepping(u,T);
u_temp = zeros(N_p,K);

% time-step choice for stable scheme
GLnodes = GaussLobattoNodes(N_fp); r_min = abs(GLnodes(2) - GLnodes(1));
dtscale = dtscale2D(); C = 0.75; dt = C * min(dtscale) * (2 / 3) * r_min / sqrt(a_advec ^ 2 + b_advec ^ 2);
% dt = 0.5 * dt; % kvuli omezeni chyby v dusledku casove diskretizace
N_steps = ceil(T / dt); dt = T / N_steps;

x_in = x(vmapI); y_in = y(vmapI);
current_time = 0;
for timestep = 1:N_steps
    for i = 1:5 % 5 stages of low-storage Runge-Kutta method
        time_RK = current_time + rk4c(i) * dt;
        
        
        % u_rhs = ODEs_RHS(u,time_RK);
        % local Lax-Friedrichs flux - strong form
        
        % computing flux on the boundaries within the elements
        flux = (a_advec * u(vmapM)) .* n_x(:) + (b_advec * u(vmapM)) .* n_y(:);
        
        % computing numerical flux
        numerical_flux = ((a_advec * (u(vmapM) + u(vmapP))) .* n_x(:) + (b_advec * (u(vmapM) + u(vmapP))) .* n_y(:) + (u(vmapM) - u(vmapP)) .* ((1 - alpha_flux) * abs(a_advec * n_x(:)+ b_advec * n_y(:)))) / 2;
        
        % imposing inflow condition
        u_inflow = BC(x_in,y_in,time_RK);
        numerical_flux(mapI) = ((a_advec * (u(vmapI) + u_inflow)) .* n_x(mapI) + (b_advec * (u(vmapI) + u_inflow)) .* n_y(mapI) + (u(vmapI) - u_inflow) .* ((1 - alpha_flux) * abs(a_advec * n_x(mapI)+ b_advec * n_y(mapI)))) / 2;
        
        % imposing outflow condition - we extrapolate u - therefore the numerical flux takes form
        numerical_flux(mapO) = (a_advec * u(vmapO)) .* n_x(mapO) + (b_advec * u(vmapO)) .* n_y(mapO);
        
        % resizing both normal and numerical fluxes to appropiate arrays
        flux = reshape(flux,N_fp * N_faces,K);
        numerical_flux = reshape(numerical_flux,N_fp * N_faces,K);
        
        % building the right-hand side of the ODE system
        [u_x,u_y] = Grad2D(u);
        u_rhs = -a_advec * u_x - b_advec * u_y + LIFT * (Fscale .* (flux - numerical_flux));
        
        
        u_temp = rk4a(i) * u_temp + dt * u_rhs; % rk4a(1) = 0, proto bude u_temp = dt * u_rhs pri kroku i = 1
        u = u + rk4b(i) * u_temp;
    end
    current_time = current_time + dt;
    
    if (rem(timestep,100) == 0)
        fprintf(2,'%i of total %i timesteps processed.\n',timestep,N_steps);
    end
end
computation_time = toc;