-
Adrien Crovato authoredAdrien Crovato authored
wSolver.cpp 9.30 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 "wMedium.h"
#include "wAssign.h"
#include "wBoundary.h"
#include "wWake.h"
#include "wLoadFunctional.h"
#include "wMshData.h"
#include "wNode.h"
#include "wElement.h"
#include "wTag.h"
#include "wResults.h"
#include "wMshExport.h"
#include <tbb/parallel_for_each.h>
#include <tbb/spin_mutex.h>
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
* @authors Adrien Crovato
*/
Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _linsol) : pbl(_pbl), linsol(_linsol), nIt(0), Cl(0), Cd(0), Cs(0), Cm(0)
{
// Say hi
std::cout << std::endl;
std::cout << "*******************************************************************************" << std::endl;
std::cout << "** \\_/ **" << std::endl;
std::cout << "** \\_______O(_)O_______/ **" << std::endl;
std::cout << "** _______________________________ **" << std::endl;
std::cout << "** ___ __ \\__ |__ __ \\__ __/ **" << std::endl;
std::cout << "** __ / / /_ /| |_ /_/ /_ / **" << std::endl;
std::cout << "** _ /_/ /_ ___ | _, _/_ / **" << std::endl;
std::cout << "** /_____/ /_/ |_/_/ |_| /_/ **" << std::endl;
std::cout << "*******************************************************************************" << std::endl;
std::cout << "** Hi! My name is DART v1.0.0-21.09 **" << std::endl;
std::cout << "** Adrien Crovato **" << std::endl;
std::cout << "** ULiege 2018-2021 **" << std::endl;
std::cout << "*******************************************************************************" << std::endl
<< std::endl;
// Default parameters
nthreads = 1;
verbose = 0;
maxIt = 50;
relTol = 1e-6;
absTol = 1e-8;
// 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
pbl->iIC->apply(phi);
for (auto dBC : pbl->dBCs)
dBC->apply(phi);
// 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;
}
Solver::~Solver()
{
// Say bye
std::cout << std::endl;
std::cout << "*******************************************************************************" << std::endl;
std::cout << "** \\_/ **" << std::endl;
std::cout << "** \\_______O(_)O_______/ **" << std::endl;
std::cout << "** ! ! **" << std::endl;
std::cout << "*******************************************************************************" << std::endl;
std::cout << std::endl;
}
/**
* @brief Dummy run
* @authors Adrien Crovato
*/
bool Solver::run()
{
throw std::runtime_error("dart::Solver::run not implemented!\n");
}
/**
* @brief Write the results
* @authors Adrien Crovato
*/
void Solver::save(MshExport *mshWriter, int n)
{
// Write files
std::cout << "Saving files... " << std::endl;
// setup 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;
// save (all mesh and boundary surface) if (n > 0)
{
mshWriter->save(pbl->msh->name + "_" + std::to_string(n), results);
for (auto sur : pbl->bnds)
sur->save(sur->groups[0]->tag->name + "_" + std::to_string(n), results);
}
else
{
mshWriter->save(pbl->msh->name, results);
for (auto sur : pbl->bnds)
sur->save(sur->groups[0]->tag->name, results);
}
}
/**
* @brief Compute total aerodynamic load coefficients on boundaries
* @authors Adrien Crovato
*/
void Solver::computeLoad()
{
// Reset coefficients and compute loads
Cl = 0;
Cd = 0;
Cs = 0;
Cm = 0;
for (auto sur : pbl->bnds)
{
// Reset
std::fill(sur->nLoads.begin(), sur->nLoads.end(), Eigen::Vector3d::Zero());
// Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element
for (auto e : sur->groups[0]->tag->elems)
{
// Build nodal pressure load
Eigen::VectorXd be = LoadFunctional::build(*e, *sur->svMap.at(e), phi, *pbl->medium->cP);
// Assembly
for (size_t i = 0; i < e->nodes.size(); ++i)
{
size_t idx = sur->nMap.at(e->nodes[i]);
sur->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)
sur->Cm = 0;
Eigen::Vector3d Cf(0, 0, 0);
for (size_t i = 0; i < sur->nodes.size(); ++i)
{
Eigen::Vector3d l = sur->nodes[i]->pos - pbl->x_ref; // lever arm
Eigen::Vector3d cf = sur->nLoads[i]; // normalized force
Cf += cf;
sur->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
sur->Cd = Cf.dot(pbl->dirD.eval());
sur->Cs = Cf.dot(pbl->dirS.eval());
sur->Cl = Cf.dot(pbl->dirL.eval());
// compute total
Cl += sur->Cl;
Cd += sur->Cd;
Cs += sur->Cs;
Cm += sur->Cm;
}
}
/**
* @brief Compute variables in the flow field
* @authors Adrien Crovato
*/void Solver::computeFlow()
{
// Compute results in volume
auto fluid = pbl->medium;
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.25)
std::cout << ANSI_COLOR_YELLOW << "Max. Mach greater than 1.25!" << ANSI_COLOR_RESET << std::endl;
}