-
Adrien Crovato authoredAdrien Crovato authored
wSolver.cpp 15.05 KiB
/*
* 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"] = φ
results.scalars_at_nodes["varPhi"] = &vPhi;
results.scalars_at_nodes["resPhi"] = &rPhi;
results.vectors_at_nodes["U"] = &U;
results.scalars_at_nodes["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;
}