solve.cpp

/* This file is a part of Numsimulation project
Developed as part of Numerical Flow Simulation course
Author: Kyrylo Sovailo */

#include <simple/simple.h>
#include <iostream>

double simple::Solver::_a(double pe) const
{
    switch (_parameters->method)
    {
        case Method::linear: return 1 - 0.5 * abs(pe);
        case Method::upwind: return 1;
        case Method::hybrid: return std::max(0.0, 1 - 0.5 * abs(pe));
        case Method::potential: pe = 1 - 0.1 * abs(pe); return std::max(0.0, pe*pe*pe*pe*pe);
        case Method::exponential: return abs(pe) / exp(abs(pe) - 1);
    }
    return 0;
}

void simple::Solver::_setup_fixed_pressure()
{
    for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
    {
        Cell *cell = *icell;
        if (cell->fixed_pressure) cell->pressure = cell->pressure_estimate = _parameters->pressure(cell->center, _time);
    }
}

void simple::Solver::_setup_fixed_velocity()
{
    for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
    {
        Face *face = *iface;
        if (face->fixed_velocity) face->velocity = face->velocity_estimate = _parameters->velocity(face->center, _time).dot(face->normal);
    }
}

void simple::Solver::_solve_velocity()
{
    double previous_residual = std::numeric_limits<double>::infinity();
    unsigned int fail_count = 0;
    unsigned int iteration = 0;
    while(true)
    {
        iteration++;
        double max_residual = -std::numeric_limits<double>::infinity();
        
        //Main calculation
        for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
        {
            Face *face = *iface;
            if (face->fixed_velocity)
            {
                //Do nothing at fixed velocity
            }
            else if (face->cells.size() < 2)
            {
                //If it's not fixed velocity, but an edge, just copy the velocity from the neighboring cell
                //Because edge cells don't have area, and don't have pressure difference
                for (unsigned int n = 0; n < face->neighbors.size(); n++)
                {
                    if (face->face_constants[n].neighbor_velocity_inwards_sign != 0)
                    {
                        (_parameters->v_gauss ? face->velocity : face->new_velocity) = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
                    }
                }
            }
            else
            {
                //Average cell, normal calculation with friction
                double a0 = _parameters->density * face->area / _parameters->step;  //a0 becomes a0_tilde till the end of the iteration
                double b = a0 * face->velocity;                                     //b becomes a0_tilde*next_velocity till the end of the iteration
                b -= face->pressure_sign * (face->cells[1]->pressure - face->cells[0]->pressure) * face->area / face->cell_cell_distance;
                for (unsigned int n = 0; n < face->neighbors.size(); n++)
                {
                    //Precalculating
                    //face velocity in the right direction, used as scalar
                    const double neighbor_velocity = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
                    //face velocity in inwards direction, used as velocity
                    const double neighbor_velocity_inwards = face->face_constants[n].neighbor_velocity_inwards_sign * neighbor_velocity;
                    const double face_face_distance = face->face_constants[n].face_face_distance;
                    const double border_length = face->face_constants[n].border_length;

                    //Kernel
                    const double d = _parameters->viscosity / face_face_distance;
                    const double pe = _parameters->density * neighbor_velocity_inwards * face_face_distance / _parameters->viscosity;
                    const double f = neighbor_velocity_inwards * _parameters->density;
                    const double a = d * border_length * _a(pe) + std::max(0.0, f * border_length);
                    a0 += a;
                    b += a * neighbor_velocity;
                }
                face->a0_tilde = a0;
                const double residual = std::abs(face->velocity * a0 - b);
                if (residual > max_residual) max_residual = residual;
                (_parameters->v_gauss ? face->velocity : face->new_velocity) = b / a0;
            }
        }

        //Applying new_velocity (Jacobi algorithm)
        if (!_parameters->v_gauss)
        {
            for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
            {
                Face *face = *iface;
                if (face->fixed_velocity) continue; //Skip fixed cells
                face->velocity = face->new_velocity;
            }
        }

        //Logics
        if (max_residual < _parameters->v_max_residual)
        {
            break;
        }
        if (iteration > _parameters->v_max_iteration)
        {
            if (_parameters->v_reaction == Reaction::error) throw std::runtime_error("simple::Solver::_solve_velocity_estimate(): Number of iterations exceeded");
            else if (_parameters->v_reaction == Reaction::warning) std::cout << "V* solver exceeded number of iterations\n";
            break;
        }
        if (max_residual > previous_residual)
        {
            if (fail_count++ > 100) break; //throw std::runtime_error("simple::Solver::_solve_velocity_estimate(): Residual does not converge");
        }
        previous_residual = max_residual;
    }
}

double simple::Solver::_solve_pressure_correction()
{
    //Precalculating c on faces
    for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
    {
        Face *face = *iface;
        if (face->fixed_velocity || face->neighbors.size() < 2) continue; //Skipped fixed velocity and edge cells faces (they will be also skipped later on)
        face->c = face->area / (face->a0_tilde * face->length);
    }

    //Precalculating e and c0_tilde on cells, setting pressure correction to zero
    double max_e = -std::numeric_limits<double>::infinity();
    Cell *max_e_cell = nullptr;
    for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
    {
        Cell *cell = *icell;
        cell->pressure_correction = 0;
        if (cell->fixed_pressure || cell->fixed_e)
        {
            //Do nothing at fixed pressure or fixed e
        }
        else
        {
            //Normal cell, calculate e and c0_tilde
            cell->e = 0;
            cell->c0_tilde = 0;
            for (unsigned int f = 0; f < cell->faces.size(); f++)
            {
                cell->e += cell->face_constants[f].face_normal_sign * cell->faces[f]->length * cell->faces[f]->velocity;
                cell->c0_tilde += cell->faces[f]->c;
            }
            if (std::abs(cell->e) > max_e) { max_e = std::abs(cell->e); max_e_cell = cell; }
        }
    }

    //Iterating
    double previous_residual = std::numeric_limits<double>::infinity();
    unsigned int fail_count = 0;
    unsigned int iteration = 0;
    while(true)
    {
        iteration++;
        double max_residual = -std::numeric_limits<double>::infinity();
        
        //Main calculation
        for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
        {
            Cell *cell = *icell;
            if (cell->fixed_pressure || cell->fixed_e)
            {
                //Do nothing at fixed pressure or fixed e
            }
            else
            {
                //Average cell, normal correction based on inflowing mass
                double sum = 0;
                for (unsigned int f = 0; f < cell->faces.size(); f++)
                {
                    if (cell->faces[f]->fixed_velocity || cell->faces[f]->neighbors.size() < 2) continue; //Skip if face is fixed velocity or edge
                    sum += cell->faces[f]->c * cell->neighbors[f]->pressure_correction;
                }
                const double residual = std::abs(cell->e + sum - cell->c0_tilde * cell->pressure_correction);
                if (residual > max_residual) max_residual = residual;
                (_parameters->p_gauss ? cell->pressure_correction : cell->new_pressure_correction) = (cell->e + sum) / cell->c0_tilde;
                if (abs(_parameters->p_gauss ? cell->pressure_correction : cell->new_pressure_correction) == std::numeric_limits<double>::infinity())
                {
                    int k = 0;
                    k++;
                }
            }
        }

        //Logics
        if (max_residual < _parameters->p_max_residual)
        {
            return max_e;
        }
        if (iteration > _parameters->p_max_iteration)
        {
            if (_parameters->p_reaction == Reaction::error) throw std::runtime_error("simple::Solver::_solve_pressure_correction(): Number of iterations exceeded");
            else if (_parameters->p_reaction == Reaction::warning) std::cout << "P' solver exceeded number of iterations\n";
            return max_e;
        }
        if (max_residual > previous_residual)
        {
            if (fail_count++ > 100) throw std::runtime_error("simple::Solver::_solve_pressure_correction(): Residual does not converge");
        }
        previous_residual = max_residual;
    }
}

void simple::Solver::_set_correct_pressure()
{
    for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
    {
        Cell *cell = *icell;
        if (cell->fixed_pressure || cell->fixed_e)
        {
            // Do nothing
        }
        else
        {
            // Normal update
            cell->pressure = cell->pressure_estimate + _parameters->beta * cell->pressure_correction;
        }
    }

    for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
    {
        Cell *cell = *icell;
        if (!cell->fixed_pressure && cell->fixed_e)
        {
            //Pressure of the cell is a sum of neighboring cells
            double sum = 0;
            double perimeter = 0;
            for (unsigned int n = 0; n < cell->neighbors.size(); n++)
            {
                if (cell->neighbors[n] != nullptr)
                {
                    sum += cell->neighbors[n]->pressure_correction * cell->faces[n]->length;
                    perimeter += cell->faces[n]->length;
                }
            }
            cell->pressure = sum / perimeter;
        }
    }
}

void simple::Solver::_set_correct_velocity()
{
    for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
    {
        Face *face = *iface;
        if (face->fixed_velocity || face->neighbors.size() < 2)
        {
            // Do nothing
        }
        else
        {
            face->velocity = face->velocity_estimate - face->pressure_sign * (face->cells[1]->pressure_correction - face->cells[0]->pressure_correction) * face->area / face->cell_cell_distance / face->a0_tilde;
        }
    }

    for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
    {
        Face *face = *iface;
        if (!face->fixed_velocity && face->neighbors.size() < 2)
        {
            // Velocity on edges is a copy of neighbor velocity
            for (unsigned int n = 0; n < face->neighbors.size(); n++)
            {
                if (face->face_constants[n].neighbor_velocity_inwards_sign != 0)
                {
                    face->velocity = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
                }
            }
        }
    }
}

void simple::Solver::_set_pressure_estimate()
{
    for (std::set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
    {
        Cell *cell = *i;
        if (cell->fixed_pressure) continue;
        cell->pressure_estimate = cell->pressure;
    }
}

void simple::Solver::_set_velocity_estimate()
{
    for (std::set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
    {
        Face *face = *i;
        if (face->fixed_velocity) continue;
        face->velocity_estimate = face->velocity;
    }
}

void simple::Solver::step()
{
    _setup_fixed_pressure();
    _setup_fixed_velocity();

    //Iteration
    double previous_max_e = std::numeric_limits<double>::infinity();
    unsigned int fail_count = 0;
    unsigned int iteration = 0;
    while(true)
    {
        iteration++;
        double max_e = _solve_pressure_correction(); //Solve P' (based on V)
        _set_correct_pressure(); //P = P* + P'
        _set_correct_velocity(); //V = V* + V'
        if (max_e < _parameters->max_e)
        {
            std::streamsize precision = std::cout.precision(4);
            std::cout << "Time step " << _time << " --> " << _time + _parameters->step << " done in " << iteration << " iterations\n";
            std::cout.precision(precision);
        }
        if (iteration > _parameters->max_iteration)
        {
            throw std::runtime_error("simple::Solver::step(): Number of iterations exceeded");
        }
        if (max_e > previous_max_e)
        {
            if (fail_count++ > 100) throw std::runtime_error("simple::Solver::step(): E does not converge");
        }
        previous_max_e = max_e;
    }

    //Solution
    _solve_velocity(); //Solve V (based on P)
    _set_pressure_estimate(); //P* = P
    _set_velocity_estimate(); //V' = V

    //Progressing
    _time += _parameters->step;
}