wSolver.cpp

/*
 * Copyright 2020 University of Liège
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "wSolver.h"
#include "wProblem.h"
#include "wFluid.h"
#include "wFreestream.h"
#include "wAssign.h"
#include "wBody.h"
#include "wWake.h"
#include "wKutta.h"
#include "wLoadFunctional.h"

#include "wMshData.h"
#include "wNode.h"
#include "wElement.h"
#include "wTag.h"
#include "wLinearSolver.h"
#include "wResults.h"
#include "wMshExport.h"

#include <tbb/parallel_for_each.h>
#include <tbb/spin_mutex.h>

#include <iomanip>
#include <fstream>

using namespace tbox;
using namespace dart;

#define ANSI_COLOR_YELLOW "\x1b[1;33m"
#define ANSI_COLOR_RESET "\x1b[0m"

/**
 * @brief Initialize the solver and perform sanity checks
 */
Solver::Solver(std::shared_ptr<Problem> _pbl,
               std::shared_ptr<LinearSolver> _linsol,
               double rTol, double aTol, int mIt,
               int nthrds, int vrb) : pbl(_pbl), linsol(_linsol),
                                      relTol(rTol), absTol(aTol), maxIt(mIt),
                                      nthreads(nthrds), verbose(vrb),
                                      nIt(0), Cl(0), Cd(0), Cs(0), Cm(0)
{
    // Say hi
    std::cout << "*******************************************************************************\n"
              << "**                                    \\_/                                    **\n"
              << "**                           \\_______O(_)O_______/                           **\n"
              << "**                      _______________________________                      **\n"
              << "**                      ___  __ \\__    |__  __ \\__  __/                      **\n"
              << "**                      __  / / /_  /| |_  /_/ /_  /                         **\n"
              << "**                      _  /_/ /_  ___ |  __ _/_  /                          **\n"
              << "**                      /_____/ /_/  |_/_/ |_| /_/                           **\n"
              << "*******************************************************************************\n"
              << "** Hi! My name is DART v1.3.0-24.12                                          **\n"
              << "** Adrien Crovato                                                            **\n"
              << "** ULiege 2018-2024                                                          **\n"
              << "*******************************************************************************\n"
              << std::endl;

    // Check the problem, pin a degree of freedom if necessary, and update element memory
    pbl->check();
    pbl->pin();
    pbl->initElems();

    // Assign a row (unknown index) for each node, and setup the local list
    rows.resize(pbl->msh->nodes.size());
    for (size_t i = 0; i < pbl->msh->nodes.size(); ++i)
        rows[i] = pbl->msh->nodes[i]->row;

    // Setup variables (unkwown vector)
    phi.resize(pbl->msh->nodes.size(), 0.);
    rPhi.resize(pbl->msh->nodes.size(), 0.);
    vPhi.resize(pbl->msh->nodes.size(), 0.);
    U.resize(pbl->msh->nodes.size(), Eigen::Vector3d::Zero());
    rho.resize(pbl->msh->nodes.size(), 0.);
    M.resize(pbl->msh->nodes.size(), 0.);
    Cp.resize(pbl->msh->nodes.size(), 0.);

    // Apply Initial and Dirichlet boundary conditions
    this->reset();

    // Kutta condition step 1: equality of normal flux through the wake
    // -> directly assemble upper wake contributions on lower wake rows, and leave the upper wake rows to zero
    // Nishida Thesis, pp. 29, Eq 2.17, Galbraith paper, pp. 5, Eq. 27
    for (auto wake : pbl->wBCs)
        for (auto p : wake->nodMap)
            rows[p.first->row] = p.second->row;

    // Display configuration
    std::cout << "--- Mesh data ---\n"
              << "Number of nodes: " << pbl->msh->nodes.size() << "\n"
              << "Number of elements: " << pbl->msh->elems.size() << "\n";
    std::cout << "--- Physical groups ---\n"
              << "Fluid on " << *pbl->fluid->tag << "\n";
    for (auto bnd : pbl->fBCs)
        std::cout << "Freestream on " << *bnd->tag << "\n";
    for (auto bnd : pbl->dBCs)
        std::cout << "Dirichlet on " << *bnd->tag << "\n";
    for (auto bnd : pbl->wBCs)
        std::cout << "Wake on " << *bnd->groups[0]->tag << "\n";
    for (auto bnd : pbl->kSCs)
        std::cout << "Kutta on " << *bnd->groups[0]->tag << "\n";
    for (auto bnd : pbl->bodies)
        std::cout << "Body on " << *bnd->groups[0]->tag << "\n";
    std::cout << "--- Freestream conditions ---\n"
              << std::setprecision(5)
              << "Angle of attack: " << pbl->alpha * 180 / 3.14159 << " deg\n"
              << "Angle of sideslip: " << pbl->beta * 180 / 3.14159 << " deg\n"
              << "Mach number: " << pbl->M_inf << "\n";
    std::cout << "--- Reference data ---\n"
              << "Reference area: " << pbl->S_ref << "\n"
              << "Reference length: " << pbl->c_ref << "\n"
              << "Reference origin: (" << pbl->x_ref(0) << ", " << pbl->x_ref(1) << ", " << pbl->x_ref(2) << ")"
              << std::endl;
}
Solver::~Solver()
{
    // Say bye
    std::cout << "\n"
              << "*******************************************************************************\n"
              << "**                                    \\_/                                    **\n"
              << "**                           \\_______O(_)O_______/                           **\n"
              << "**                                  !     !                                  **\n"
              << "*******************************************************************************\n"
              << std::endl;
}

/**
 * @brief Set or reset initial conditions
 */
void Solver::reset()
{
    // ICs
    pbl->iIC->apply(phi);
    // BCs
    for (auto dBC : pbl->dBCs)
        dBC->apply(phi);
}

/**
 * @brief Dummy run
 */
Status Solver::run()
{
    throw std::runtime_error("dart::Solver::run not implemented!\n");
}

/**
 * @brief Write the results
 */
void Solver::save(MshExport *mshWriter, std::string const &suffix)
{
    // Write files
    std::cout << "Saving files "
              << std::fixed << std::setprecision(2)
              << "(Mach " << pbl->M_inf << ", "
              << pbl->alpha * 180 / 3.14159 << "deg AoA, "
              << pbl->beta * 180 / 3.14159 << "deg AoS)"
              << std::endl;
    // save results
    Results results;
    results.scalars_at_nodes["phi"] = &phi;
    results.scalars_at_nodes["varPhi"] = &vPhi;
    results.scalars_at_nodes["resPhi"] = &rPhi;
    results.vectors_at_nodes["U"] = &U;
    results.scalars_at_nodes["rho"] = &rho;
    results.scalars_at_nodes["Mach"] = &M;
    results.scalars_at_nodes["Cp"] = &Cp;
    mshWriter->save(pbl->msh->name + suffix, results);
    // save aerodynamic loads
    std::cout << "writing file: aeroloads" << suffix << ".dat... " << std::flush;
    std::ofstream outfile;
    outfile.open("aeroloads" + suffix + ".dat");
    // total
    outfile << "# Total - " << pbl->msh->name << std::endl;
    outfile << std::fixed << std::setprecision(12)
            << "Cl = " << Cl << "\n"
            << "Cd = " << Cd << "\n"
            << "Cs = " << Cs << "\n"
            << "Cm = " << Cm << "\n";
    // breakdown
    for (auto b : pbl->bodies)
    {
        outfile << "# Body - " << b->groups[0]->tag->name << " (" << b->groups[0]->tag->elems.size() << " elements)" << std::endl;
        outfile << std::fixed << std::setprecision(12)
                << "Cl = " << b->Cl << "\n"
                << "Cd = " << b->Cd << "\n"
                << "Cs = " << b->Cs << "\n"
                << "Cm = " << b->Cm << "\n";
    }
    outfile.close();
    std::cout << "done." << std::endl;
}

/**
 * @brief Compute refinment sensor and write new mesh sizes
 * @todo Remesh internally
 * @todo Investigate combination of gradient sensor with Laplacian sensor
 * @todo Investigate effect of user-defined parameters (thresholds, refinement/coarsening factors, termination, etc.)
 * @todo Add capability to define region of interest (LE, TE, exclude tip, etc.)
 */
bool Solver::refine(double minSize, tbox::MshExport *mshWriter, std::string const &suffix)
{
    // Compute deltaRho at cell in the streamwise direction
    std::vector<double> dRhoE(pbl->msh->elems.size());
    std::vector<double> dRho(pbl->msh->nodes.size()), mSize(pbl->msh->nodes.size());
    tbb::parallel_for_each(pbl->fluid->adjMap.begin(), pbl->fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> &p) {
        dRhoE[p.first->no - 1] = std::abs(pbl->fluid->rho->eval(*(p.first), phi, 0) - pbl->fluid->rho->eval(*(p.second[0]), phi, 0));
    });
    // Average deltaRho and compute mesh size at nodes
    tbb::parallel_for_each(pbl->fluid->neMap.begin(), pbl->fluid->neMap.end(), [&](std::pair<Node *, std::vector<Element *>> p) {
        dRho[p.first->row] = 0.;
        mSize[p.first->row] = 0.;
        for (auto e : p.second)
        {
            dRho[p.first->row] += dRhoE[e->no - 1];
            if (pbl->nDim == 2)
                mSize[p.first->row] += std::sqrt(2 / std::sqrt(3) * e->getDetJ(0));
            else
                mSize[p.first->row] += std::cbrt(2 / std::sqrt(2) * e->getDetJ(0));
        }
        dRho[p.first->row] /= p.second.size();
        mSize[p.first->row] /= p.second.size();
    });

    // Sort the mesh size of the element according the their deltaRho
    std::vector<std::pair<int, double>> sorter(pbl->msh->nodes.size());
    tbb::parallel_for_each(pbl->msh->nodes.begin(), pbl->msh->nodes.end(), [&](Node *n) {
        sorter[n->row] = std::pair<int, double>(n->row, dRho[n->row]);
    });
    std::sort(sorter.begin(), sorter.end(), [&](const std::pair<int, double> &a, const std::pair<int, double> &b) {
        return a.second > b.second;
    });
    // Coarsen the 30% elements which have lowest dRho, and refine the 30% that have highest dRho
    size_t refBnd = static_cast<size_t>(std::floor(0.3 * sorter.size()));
    size_t crsBnd = static_cast<size_t>(std::floor(0.7 * sorter.size()));
    tbb::parallel_for(tbb::blocked_range<size_t>(0, sorter.size()), [&](tbb::blocked_range<size_t> rng) {
        for (size_t i = rng.begin(); i < rng.end(); ++i)
        {
            // Manually prevent refine in certain regions
            // Node *n = pbl->msh->nodes[sorter[i].first];
            // if (n->pos(1) > 0.98 * 14 && n->pos(1) < 1.02 * 14)
            //    if (mSize[sorter[i].first] < 0.045)
            //        continue;
            if (i < refBnd)
                mSize[sorter[i].first] *= 0.5;
            else if (i > crsBnd)
                mSize[sorter[i].first] *= 2.0;
            else
                mSize[sorter[i].first] *= 1.0;
        }
    });

    // Write sensor and new mesh sizes
    Results rSens;
    rSens.scalars_at_nodes["sensor"] = &dRho;
    mshWriter->saveList(pbl->msh->name + "_sensor" + suffix, rSens);
    Results rSize;
    rSize.scalars_at_nodes["size"] = &mSize;
    mshWriter->saveList(pbl->msh->name + "_size" + suffix, rSize);

    // Stop the refinement loop if the mesh size is lower that the specified target
    std::cout << std::fixed << std::setprecision(5) << "min. size: " << *std::min_element(mSize.begin(), mSize.end()) << std::endl;
    if (*std::min_element(mSize.begin(), mSize.end()) < minSize)
        return true;
    else
        return false;
}

/**
 * @brief Compute total aerodynamic load coefficients on boundaries
 */
void Solver::computeLoad()
{
    // Reset coefficients and compute loads
    Cl = 0;
    Cd = 0;
    Cs = 0;
    Cm = 0;
    for (auto bnd : pbl->bodies)
    {
        // Reset
        std::fill(bnd->nLoads.begin(), bnd->nLoads.end(), Eigen::Vector3d::Zero());

        // Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element
        for (auto e : bnd->groups[0]->tag->elems)
        {
            // Build nodal pressure load
            Eigen::VectorXd be = LoadFunctional::build(*e, *bnd->svMap.at(e), phi, *pbl->fluid->cP);

            // Assembly
            for (size_t i = 0; i < e->nodes.size(); ++i)
            {
                size_t idx = bnd->nMap.at(e->nodes[i]);
                bnd->nLoads[idx] += be.block<3, 1>(i * 3, 0);
            }
        }
        // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
        // force coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
        bnd->Cm = 0;
        Eigen::Vector3d Cf(0, 0, 0);
        for (size_t i = 0; i < bnd->nodes.size(); ++i)
        {
            Eigen::Vector3d l = bnd->nodes[i]->pos - pbl->x_ref; // lever arm
            Eigen::Vector3d cf = bnd->nLoads[i];                 // normalized force
            Cf += cf;
            bnd->Cm += (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
        }
        // rotate to flow direction
        bnd->Cd = Cf.dot(pbl->dirD.eval());
        bnd->Cs = Cf.dot(pbl->dirS.eval());
        bnd->Cl = Cf.dot(pbl->dirL.eval());
        // compute total
        Cl += bnd->Cl;
        Cd += bnd->Cd;
        Cs += bnd->Cs;
        Cm += bnd->Cm;
    }
}

/**
 * @brief Compute variables in the flow field
 */
void Solver::computeFlow()
{
    // Compute results in volume
    auto fluid = pbl->fluid;
    tbb::parallel_for_each(fluid->neMap.begin(), fluid->neMap.end(), [&](std::pair<Node *, std::vector<Element *>> p) {
        Node *n = p.first;
        // Perturbation potential
        vPhi[n->row] = phi[n->row] - fluid->phiInf->eval(n->pos);
        // Velocity (averaged at nodes)
        U[n->row] = Eigen::Vector3d::Zero();
        for (auto e : p.second)
        {
            Eigen::VectorXd grad = e->computeGradient(phi, 0);
            for (int i = 0; i < grad.size(); ++i)
                U[n->row](i) += grad(i);
        }
        U[n->row] /= static_cast<double>(p.second.size());
        // Density (averaged at nodes)
        rho[n->row] = 0.;
        for (auto e : p.second)
            rho[n->row] += fluid->rho->eval(*e, phi, 0);
        rho[n->row] /= p.second.size();
        // Mach number (averaged at nodes)
        M[n->row] = 0.;
        for (auto e : p.second)
            M[n->row] += fluid->mach->eval(*e, phi, 0);
        M[n->row] /= p.second.size();
        // Pressure coefficient (averaged at nodes)
        Cp[n->row] = 0.;
        for (auto e : p.second)
            Cp[n->row] += fluid->cP->eval(*e, phi, 0);
        Cp[n->row] /= p.second.size();
    });

    // Check maximum Mach number
    if (*std::max_element(M.begin(), M.end()) >= 1.3 && verbose > 0)
        std::cout << ANSI_COLOR_YELLOW << "Max. Mach greater than 1.3." << ANSI_COLOR_RESET << std::endl;
}