clc; clear; close all;

Globals1D;

tic;
Problem;

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

StartUp1D; % computing local operators and mesh information

[mapI, vmapI, mapO, vmapO] = BuildBCMaps();

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

u = IC(x); % evaluating initial condition
preprocessing_time = toc;

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

% time-step choice for stable scheme
CLF = 0.75; dt = (CLF / abs(a_advec)) * min(abs(x(1,:) - x(2,:)));
% dt = 0.25 * dt; % kvuli omezeni chyby v dusledku casove diskretizace
N_steps = ceil(T / dt); dt = T / N_steps;

x_in = x(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(:);
        
        % computing numerical flux
        numerical_flux = ((a_advec * (u(vmapM) + u(vmapP))) .* n_x(:) + (u(vmapM) - u(vmapP)) .* ((1 - alpha_flux) * abs(a_advec * n_x(:)))) / 2;
        
        % imposing inflow condition
        u_inflow = BC(x_in,time_RK);
        numerical_flux(mapI) = ((a_advec * (u(vmapI) + u_inflow)) .* n_x(mapI) + (u(vmapI) - u_inflow) .* ((1 - alpha_flux) * abs(a_advec * n_x(mapI)))) / 2;
        
        % imposing outflow condition - we extrapolate u - therefore the numerical flux takes form
        numerical_flux(mapO) = (a_advec * u(vmapO)) .* n_x(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_rhs = -a_advec * r_x .* (D_r * u) + 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;