//#####################################################################
// Copyright 2002-2004, Ronald Fedkiw.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Advection/ADVECTION_CONSERVATIVE_WENO.h>
#include <Core/Advection/ADVECTION_SEPARABLE_UNIFORM.h>
#include <Core/Math_Tools/minmax.h>
namespace PhysBAM{
//#####################################################################
// Function Advection_Solver
//#####################################################################
// finds (uZ)_x with Local Lax Friedrichs
template<class T,class T2> void ADVECTION_CONSERVATIVE_WENO<T,T2>::
Advection_Solver(const int m,const T dx,const ARRAY<T2,IV>& Z,const ARRAY<T,IV>& u,ARRAY<T2,IV>& u_Zx,
    ARRAY<T2>* flux_input)
{
    ARRAY<T2,IV> DZ0(RANGE<IV>(IV(-3),IV(m+3))); // 1st divided difference
    for(int i=-3;i<m+3;i++) DZ0(i)=Z(i);
    ARRAY<T2,IV> DUZ0(RANGE<IV>(IV(-3),IV(m+3))); // 1st divided difference
    for(int i=-3;i<m+3;i++) DUZ0(i)=u(i)*Z(i);

    if(compute_epsilon){
        epsilon=1e-6*sqr(DUZ0.Max_Abs()); // only DUZ used to find epsilon 
        if(epsilon == 0) epsilon=1e-6;} // epsilon=0 implies all v_i=0 and all u_Zx=0

    ARRAY<T2> flux_raw;
    if(!flux_input) flux_raw.Resize(m+1); // flux is to the left of each point
    ARRAY<T2>& flux=flux_input?*flux_input:flux_raw;

    for(int i=-1;i<m;i++){
        T2 alpha=maxabs(u(i),u(i+1));           
        T2 fluxleft=WENO(DUZ0(i-2)+alpha*DZ0(i-2),DUZ0(i-1)+alpha*DZ0(i-1),DUZ0(i)+alpha*DZ0(i),DUZ0(i+1)+alpha*DZ0(i+1),DUZ0(i+2)+alpha*DZ0(i+2),epsilon);
        T2 fluxright=WENO(DUZ0(i+3)-alpha*DZ0(i+3),DUZ0(i+2)-alpha*DZ0(i+2),DUZ0(i+1)-alpha*DZ0(i+1),DUZ0(i)-alpha*DZ0(i),DUZ0(i-1)-alpha*DZ0(i-1),epsilon);
        flux(i+1)=.5*(fluxleft+fluxright);}

    T one_over_dx=1/dx;

    for(int i=0;i<m;i++) u_Zx(i)=(flux(i+1)-flux(i))*one_over_dx;
}

template class ADVECTION_CONSERVATIVE_WENO<VECTOR<double,1>,double>;
template class ADVECTION_CONSERVATIVE_WENO<VECTOR<double,2>,double>;
template class ADVECTION_CONSERVATIVE_WENO<VECTOR<double,3>,double>;
template class ADVECTION_CONSERVATIVE_WENO<VECTOR<float,2>,float>;
template class ADVECTION_CONSERVATIVE_WENO<VECTOR<float,3>,float>;

}
