//-----------------------------------------------------------------------------
// NastSolverPoisson2dSOR.cpp
//-----------------------------------------------------------------------------
//
//  Copyright (C) 1998 Technische Universitaet Muenchen, Germany
//  written by Bernhard Brueck
//
//  This file is part of Nast++
//
//-----------------------------------------------------------------------------
//  CNastSolverPoisson2dSOR ist eine Implementierung des Poissonloesers unter
//  der Verwendung des SOR-Verfahrens.
//-----------------------------------------------------------------------------
//  Changes:
//

#include "NastSolverPoisson2dSOR.h"
#include "NastDebug.h"
#include "NastStaggeredGrid2d.h"
#include "NastCell2d.h"
#include "NastArray.h"		// fuer NAST_DUMP_VISUAL_INTERVALL

//-----------------------------------------------------------------------------
//                    Konstruktor + Destruktor
//-----------------------------------------------------------------------------

CNastSolverPoisson2dSOR::CNastSolverPoisson2dSOR( int    itermax,  
        					  double eps,   
        					  double omega,
        					  bool   bUseRedBlack,
        					  int    nDebug )
    :m_rdx2(0.0),		// wird in solve festgelegt
     m_rdy2(0.0),               // wird in solve festgelegt
     m_omega(omega),		// Relaxationsparameter fuer SOR
     m_pGrid(0),		// Zeiger auf das Gitter
     m_pRhs(0),			// Zeiger auf die rechte Seite
     m_itermax(itermax),	// max Anzahl von Iterationen
     m_eps(eps),	       	// fuer die Abbruchbedingung
     m_bUseRedBlack( bUseRedBlack ),
     m_nDebug( nDebug )
{
    // Test ob die Parameter in vernuenftigen Bereichen liegen
    NAST_ASSERT( itermax > 0 );	
    NAST_ASSERT( eps > 0 );	
    NAST_ASSERT( omega >  0 );
    NAST_ASSERT( omega <= 2 );
}


//  Destruktor
CNastSolverPoisson2dSOR::~CNastSolverPoisson2dSOR()
{
    NAST_ASSERT_VALID( this );
    // vorerst nichts zu tun
}


//-------------------------------------------------------------------------
//                            Hilfsfunktionen
//-------------------------------------------------------------------------
    

inline CNastGrid2d& CNastSolverPoisson2dSOR::rhs() const
{
    NAST_ASSERT( m_pRhs );
    return *m_pRhs;
}

inline CNastStaggeredGrid2d& CNastSolverPoisson2dSOR::grid() const
{
    NAST_ASSERT( m_pGrid );
    return *m_pGrid;
}


// Residium fuer die Zelle i,j berechnen
inline double 
CNastSolverPoisson2dSOR::residium( int i, int j) const
{
    NAST_ASSERT( grid().cell(i,j).isFluid());

    double diff_x = 0;
    double diff_y = 0;

    if( grid().cell(i,j).isInsideFluid() ) // Fluidzelle umgeben mit Fluidzellen ?
    {   
        const double tmp =  -2 * grid().P(i,j);
        diff_x = (grid().P(i-1,j  ) + tmp + grid().P(i+1,j  ));
        diff_y = (grid().P(i  ,j-1) + tmp + grid().P(i  ,j+1));
    }
    else
    {
        double p = grid().P(i,j);
        if( grid().cell(i,j).isFluid( CNastCell2d::SOUTH ))  diff_y += grid().P(i  ,j-1) - p;
        if( grid().cell(i,j).isFluid( CNastCell2d::NORTH ))  diff_y += grid().P(i  ,j+1) - p;
        if( grid().cell(i,j).isFluid( CNastCell2d::WEST  ))  diff_x += grid().P(i-1,j  ) - p;
        if( grid().cell(i,j).isFluid( CNastCell2d::EAST  ))  diff_x += grid().P(i+1,j  ) - p;
    }

    return diff_x * m_rdx2 + diff_y * m_rdy2 - rhs()(i,j);
}


// zur Fehlersuche Residium anzeigen
void
CNastSolverPoisson2dSOR::show_residium() const
{
    CNastGrid2d res( rhs().box(),
        	     rhs().sizeI(), 
        	     rhs().sizeJ());
    res.setAll(0.0);

    for( int i = 1; i <= grid().sizeI(); i++)
        for( int j = 1; j <= grid().sizeJ(); j++)
            if( grid().cell(i,j).isFluid())
        	res(i,j) = residium(i,j);

    NAST_DUMP_VISUAL( res );
    NAST_DUMP_VISUAL( grid().gridP() );
    NAST_BREAK();
}


//  das Residium  ueber eine ganze Zeile
double 
CNastSolverPoisson2dSOR::residium( int i ) const
{
    double sum = 0;		// Summe der Quadrate

    for( int j = 1; j <= grid().sizeJ(); j++)
        if( grid().cell(i,j).isFluid())
        {
            double r = residium( i, j );
            sum += r*r;
        }

    return sum;
}

//  das Residium fuer die Abbruchbedingung berechnen
double 
CNastSolverPoisson2dSOR::residium() const
{
    double sum = 0;		// Summe der Quadrate

    for( int i = 1; i <= grid().sizeI(); i++)
        sum += residium(i);
    
    return sqrt(sum / grid().countFluidCells());
}


// ein Relaxationschritt ueber eine komplette Zeile
double
CNastSolverPoisson2dSOR::relax_red_black( int i, bool toggleRedBlack ) const
{
    double fac = m_omega /(2.0*(m_rdx2+m_rdy2));
    double sum = 0 ;

    int nAdd = 0;
    if( toggleRedBlack ) 
        nAdd = 1;

    for( int j = 1 + (i + nAdd) % 2; j <= grid().sizeJ(); j += 2)
    {
        if( grid().cell(i,j).isFluid() )
        {
            double res = residium( i, j );
            sum += res * res;
            grid().setP( i, j , grid().P(i,j) + res * fac );
        }
    }

    return sum;
}


//  Relaxationschritt ueber alle Zellen
double 
CNastSolverPoisson2dSOR::relax_red_black() const
{
    double res = 0;

    res += relax_red_black( 1, true );
    for( int i = 1; i < grid().sizeI(); i++)
    {
        res += relax_red_black( i+1, true);
        res += relax_red_black( i,   false);
    }
    res += relax_red_black(  grid().sizeI(), false);

    return sqrt(res / grid().countFluidCells());
}


void 
CNastSolverPoisson2dSOR::relax( int i ) const
{
    double fac = m_omega / ( 2.0 * (m_rdx2 + m_rdy2) );

    for( int j = 1; j <= grid().sizeJ(); j++ )
    {
        if( grid().cell(i,j).isFluid() )
        {
            double add = fac * residium( i, j );
            grid().setP( i, j , grid().P(i,j) + add );
        }
    }
}

// Hier werden zwei Relaxationsschritte ausgefuehrt und danach das
// Residium berechnet. 
double
CNastSolverPoisson2dSOR::relax() const
{
    double sum = 0;

    relax( 1 );
    for( int i = 1; i <= grid().sizeI() - 1; i++)
    {
        relax( i+1 );
        relax( i );
        if( i > 1 ) sum += residium(i-1);
    }
    relax( grid().sizeI() );

    sum += residium( grid().sizeI() - 1 );
    sum += residium( grid().sizeI()     );

    return sum;
}

// Schleife wird solange ausgefuehrt bis das Residium klein genug ist
void 
CNastSolverPoisson2dSOR::loop() 
{
    double res = 0;
    int    n = 0;

    if( m_nDebug == 2 ) 
        NAST_DUMP_VISUAL( rhs() );
    do
    {
        if( m_bUseRedBlack )
        {
            res = relax_red_black();
            n++;
        }
        else
        {
            res = relax();
            n += 2;
        }
        
        // Debugausgabe gewuenscht ?
        if( m_nDebug == 2 ) 
        {
            show_residium();
            NAST_DUMP_VISUAL_INTERVALL( res, 20);
        }
    }
    while( res > m_eps &&	// solange bis das Residium klein genug ist,
           n < m_itermax );	// oder die max. Anzahl der Schritte erreicht wurde

    CNastDebug::dumpContextLogfile() << "\tn:" << n;
    CNastDebug::dumpContextLogfile() << "\tres:" << res << "\n";
        
    if( m_nDebug == 1 ) show_residium();
}

//-------------------------------------------------------------------------
//                            Memberfunktionen
//-------------------------------------------------------------------------

void 
CNastSolverPoisson2dSOR::solve( CNastStaggeredGrid2d &gridP,
        			CNastGrid2d          &gridRhs ) 
{
    NAST_ASSERT_VALID( this );

    // Gitter und rechte Seite merken
    m_pGrid = &gridP;
    m_pRhs  = &gridRhs;

    //  Variablen belegen
    NAST_ASSERT( grid().dx() > 0 );
    NAST_ASSERT( grid().dy() > 0 );

    m_rdx2  = 1 / (grid().dx() * grid().dx());
    m_rdy2  = 1 / (grid().dy() * grid().dy());

    //  Summe ueber die gesamte rechte Seite muss 0 ergeben
#   ifdef NAST_DEBUG
    double sum_neg = 0;
    double sum_pos = 0;

    for( int i = 0; i < rhs().sizeI(); i++)
        for( int j = 0; j < rhs().sizeJ(); j++)
        {
            double r = rhs()(i,j);
            if( r > 0 )
        	sum_pos += r;
            else
        	sum_neg += r;
        }

    // um die Probleme mit Ungenauigkeiten zu vermeiden wird
    // nicht getestet ob die summe gleich null ist sondern
    // ob sie mehrere Groessenordnungen unter der Summe aller
    // positven Elemente ist
    NAST_ASSERT( sum_pos + sum_neg < sum_pos * 1E-12 );
#   endif 

    // und rein in die Hauptschleife
    loop();

    // und danach wieder aufraeumen
    m_rdx2  = 0;
    m_rdy2  = 0;
    m_pGrid = 0;
    m_pRhs  = 0;
}    


// ----------------------------------------------------------------------------
//                                    Debug
// ----------------------------------------------------------------------------
void 
CNastSolverPoisson2dSOR::debugDump( CNastDumpContext &dumpContext ) const
{
    CNastSolverPoisson2d::debugDump( dumpContext );
    dumpContext << "CNastSolverPoisson2dSOR" << "\n";
    dumpContext << "\t" << "m_eps     = " << m_eps     << "\n";
    dumpContext << "\t" << "m_omega   = " << m_omega   << "\n";
    dumpContext << "\t" << "m_itermax = " << m_itermax << "\n";
    dumpContext << "\t" << "\n";
    dumpContext << "\t" << "m_rdx2    = " << m_rdx2    << "\n";
    dumpContext << "\t" << "m_rdy2    = " << m_rdy2    << "\n";
}


void 
CNastSolverPoisson2dSOR::assertValid() const
{
    // zuerst einmal AssertValid der Basisklasse aufrufen
    CNastSolverPoisson2d::assertValid(); 
    NAST_ASSERT( m_eps     >  0 );
    NAST_ASSERT( m_itermax >  0 );
    NAST_ASSERT( m_omega   >  0 );
    NAST_ASSERT( m_omega   <= 2 );
}