Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • Paul.Dechamps/dartflo
1 result
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 403 additions and 253 deletions
dart/src/wKuttaElement.cpp
View file @ 5ba134b1
......@@ -26,28 +26,59 @@ KuttaElement::KuttaElement(size_t _no, Element *&_surE, Element *&_volE,
{
// Sanity checks
if (surE->type() == ELTYPE::LINE2)
if (surE->type() == ElType::LINE2)
nDim = 2;
else if (surE->type() == ElType::TRI3)
nDim = 3;
else
{
std::stringstream err;
err << "KuttaElement only implemented for surface elements of type "
<< ELTYPE::LINE2 << " (" << surE->type() << " was given)!\n";
<< ElType::LINE2 << " or " << ElType::TRI3
<< " (" << surE->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
if (volE->type() == ELTYPE::TRI3)
if (volE->type() == ElType::TRI3 || volE->type() == ElType::TETRA4)
nCol = volE->nodes.size();
else
{
std::stringstream err;
err << "KuttaElement only implemented for volume elements of type "
<< ELTYPE::TRI3 << " (" << volE->type() << " was given)!\n";
<< ElType::TRI3 << " or " << ElType::TETRA4
<< " (" << volE->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
// Count TE nodes
mapNodes();
nRow = teN.size();
}
/**
* @brief Map the wake nodes to their local volume node indices
*/
void KuttaElement::mapNodes()
{
// Find local volume index of contributing surface nodes
for (auto p : teN)
{
for (size_t i = 0; i < volE->nodes.size(); ++i)
{
if (surE->nodes[p.first] == volE->nodes[i])
{
vsMap[p.first] = i;
break;
}
}
}
// Check that all contributing surface nodes have been mapped
if (vsMap.size() != teN.size())
{
std::stringstream err;
err << "dart::KuttaElement: could not map some contributing nodes from surface element " << *surE << " to " << *volE << " (number of mapped nodes: " << vsMap.size() << ")!\n";
throw std::runtime_error(err.str());
}
}
void KuttaElement::write(std::ostream &out) const
{
out << "KuttaElement #" << no
......
dart/src/wKuttaElement.h
View file @ 5ba134b1
......@@ -21,6 +21,7 @@
#include "wObject.h"
#include "wElement.h"
#include <vector>
#include <map>
namespace dart
{
......@@ -29,16 +30,12 @@ namespace dart
* @brief Super-element to apply Kutta condition
*
* Made up of one surface element at the trailing edge and the connected volume
* element. Only Line2 (surface) and Tri3 (volume) elements are curently
* implemented.
* Since Tri3 are used, the evaluation of the volume integrals is performed at
* dummy Gauss point #0. If volume elements with non-constant shape function
* derivatives are used, those derivates should be recomputed at the true
* surface element Gauss points.
* Kutta works well in 2D, but fails in 3D unless the geometry follows strict
* restrictions (thin and sharp TE, simple planform, same mesh elements on TE
* suction and pressure side, ...). Kutta contribution is therefore disabled for
* 3D cases.
* element. Only Line2/Tri3 (surface) and Tri3/Tetra4 (volume) elements are
* currently implemented.
* Since Tri3 or Tetra4 are used, the evalutaion of the volume integrals is
* performed at dummy Gauss point #0. If volume elements with non-constant shape
* function derivatives are used, those derivates should be recomputed at the
* true surface element Gauss points.
* @authors Adrien Crovato
*/
class DART_API KuttaElement : public fwk::wObject
......@@ -54,14 +51,18 @@ public:
size_t nDim; ///< dimension of volume element
std::vector<std::pair<size_t, tbox::Node *>> teN; ///< Map of local node indices and trailing edge nodes
std::map<size_t, size_t> vsMap; ///< map between surface node and volume node indices
int sign; ///< +1 if upper element, -1 if lower element
KuttaElement(size_t _no, tbox::Element *&_surE, tbox::Element *&_volE, std::vector<std::pair<size_t, tbox::Node *>> &_teN, int _sign);
virtual void write(std::ostream &out) const override;
#endif
private:
void mapNodes();
};
} // namespace dart
#endif //WKUTTAELEMENT_H
#endif // WKUTTAELEMENT_H
dart/src/wKuttaResidual.cpp
View file @ 5ba134b1
......@@ -27,18 +27,33 @@ using namespace dart;
/**
* @brief Build the residual matrix for a fixed-point iteration, on one Kutta element
*
* A = \int psi * v_u*grad_phi dS
* A = 1/S * \int (psi + grad_psi*v_inf) * v*grad_phi dS
* NB: v_inf = [1, 0, 0]
*/
Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<double> const &phi)
{
// Get shape functions, Jacobian and Gauss points
// Get shape functions and Gauss points
Cache &sCache = ke.surE->getVCache();
Gauss &sGauss = sCache.getVGauss();
Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0);
// Get Jacobian
Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
// Get size
size_t nRow = ke.nRow;
size_t nDim = ke.nDim;
// Stabilization (disabled in 2D)
Eigen::VectorXd dSfp = Eigen::VectorXd::Zero(nRow);
if (nDim == 3)
{
Eigen::MatrixXd dSf(nRow, nDim);
for (size_t i = 0; i < nRow; ++i)
for (size_t j = 0; j < nDim; ++j)
dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first));
dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
}
// Velocity
Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getJinv(0) * ke.volE->getVCache().getDsf(0);
Eigen::RowVectorXd Aup = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).transpose() * iJ * dSfV;
// Build
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(nRow, ke.nCol);
......@@ -48,7 +63,7 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
sfp(i) = sf(ke.teN[i].first);
sfp(i) = sf(ke.teN[i].first) + dSfp(i);
// wake contribution
A += ke.sign * sfp * Aup * sGauss.getW(k) * ke.surE->getDetJ(k);
}
......@@ -58,19 +73,33 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
/**
* @brief Build the residual vector, on one Kutta element
*
* b = \int psi * (grad_phi)^2 dS
* b = 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS
* NB: v_inf = [1, 0, 0]
*/
Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double> const &phi)
{
// Get shape functions, Jacobian and Gauss points
// Get shape functions and Gauss points
Cache &sCache = ke.surE->getVCache();
Gauss &sGauss = sCache.getVGauss();
Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0);
// Get Jacobian
Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
// Get size
size_t nRow = ke.nRow;
size_t nDim = ke.nDim;
// Stabilization (disabled in 2D)
Eigen::VectorXd dSfp = Eigen::VectorXd::Zero(nRow);
if (nDim == 3)
{
Eigen::MatrixXd dSf(nRow, nDim);
for (size_t i = 0; i < nRow; ++i)
for (size_t j = 0; j < nDim; ++j)
dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first));
dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
}
// Velocity
Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0);
double dPhi2 = dPhi.squaredNorm();
double dPhi2 = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).squaredNorm();
// Build
Eigen::VectorXd b = Eigen::VectorXd::Zero(nRow);
......@@ -80,7 +109,7 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
sfp(i) = sf(ke.teN[i].first);
sfp(i) = sf(ke.teN[i].first) + dSfp(i);
// wake contribution
b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * ke.surE->getDetJ(k);
}
......@@ -90,7 +119,8 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
/**
* @brief Build the gradient of the residual vector with respect to the flow variable (jacobian), on one Kutta element
*
* A = \int psi * 2*(grad_phi_u)*dgrad_phi_u dS
* A = 1/S * \int (psi + grad_psi*v_inf) * 2*(grad_phi_u)*dgrad_phi_u dS
* NB: v_inf = [1, 0, 0]
*/
Eigen::MatrixXd KuttaResidual::buildGradientFlow(KuttaElement const &ke, std::vector<double> const &phi)
{
......@@ -102,26 +132,67 @@ Eigen::MatrixXd KuttaResidual::buildGradientFlow(KuttaElement const &ke, std::ve
/**
* @brief Build the gradient of the residual vector with respect to the nodes, on one Kutta element
*
* A = d( \int psi * (grad_phi)^2 dS )
* = \int psi * d(grad_phi)^2 dS
* + \int psi * (grad_phi)^2 ddS
* A = d( 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS )
* = d1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS
* + 1/S * \int dgrad_psi*v_inf * (grad_phi)^2 dS
* + 1/S * \int (psi + grad_psi*v_inf) * d(grad_phi)^2 dS
* + 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 ddS
* NB: v_inf = [1, 0, 0]
*/
std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(KuttaElement const &ke, std::vector<double> const &phi)
{
// Get shape functions and Gauss points
Cache &surCache = ke.surE->getVCache();
Gauss &surGauss = surCache.getVGauss();
Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0);
// Get Jacobian
Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
std::vector<Eigen::MatrixXd> const &dJ = ke.volE->getGradJ(0);
// Get surface and gradient
double surf = ke.surE->getVol();
std::vector<double> const &dSurf = ke.surE->getGradVol();
// Get size
size_t nRow = ke.nRow;
size_t nDim = ke.nDim;
size_t nCol = ke.nCol;
size_t nSur = ke.surE->nodes.size();
// velocity
// Stabilization (disabled in 2D)
Eigen::VectorXd gradSfp = Eigen::VectorXd::Zero(nRow);
std::vector<Eigen::VectorXd> dGradSfV(nDim * nCol, Eigen::VectorXd::Zero(nRow));
std::vector<Eigen::VectorXd> dGradSfS(nDim * nSur, Eigen::VectorXd::Zero(nRow));
if (nDim == 3)
{
// term
Eigen::VectorXd vi = Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
double sqSurf = sqrt(surf);
Eigen::MatrixXd dSf(nRow, nDim);
for (size_t i = 0; i < nRow; ++i)
for (size_t m = 0; m < nDim; ++m)
dSf(i, m) = dSfV(m, ke.vsMap.at(ke.teN[i].first));
gradSfp = 0.5 * sqSurf * dSf * iJ.transpose() * vi;
// volume gradient
for (size_t i = 0; i < nCol; ++i)
{
for (size_t m = 0; m < nDim; ++m)
{
size_t idx = i * nDim + m;
dGradSfV[idx] = 0.5 * sqSurf * dSf * (-iJ * dJ[idx] * iJ).transpose() * vi; // S * dgrad_psi
}
}
// surface gradient
for (size_t i = 0; i < nSur; ++i)
{
for (size_t m = 0; m < nDim; ++m)
{
size_t idx = i * nDim + m;
dGradSfS[idx] = 0.5 * 0.5 / sqSurf * dSurf[idx] * dSf * iJ.transpose() * vi; // dS * grad_psi
}
}
}
// Velocity and normalization
Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0);
double nrm = 1 / surf;
// Build
Eigen::MatrixXd Av = Eigen::MatrixXd::Zero(nRow, nDim * nCol);
......@@ -136,7 +207,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
Eigen::VectorXd const &sf = surCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
sfp(i) = sf(ke.teN[i].first);
sfp(i) = sf(ke.teN[i].first) + gradSfp(i);
// Volume
for (size_t i = 0; i < nCol; ++i)
......@@ -145,7 +216,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
{
size_t idx = i * nDim + m;
// gradient of velocity jump
Av.col(idx) += ke.sign * w * sfp * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // psi * dgrad_phi^2
Av.col(idx) += ke.sign * nrm * w * sfp * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // 1/S * (psi + grad_psi) * dgrad_phi^2 * detJ
// gradient of stabilization term (volume)
Av.col(idx) += ke.sign * nrm * w * dGradSfV[idx] * dPhi.dot(dPhi) * detJ; // 1/S * dgrad_psi * grad_phi^2 * detJ
}
}
// Surface
......@@ -154,8 +227,12 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
for (size_t m = 0; m < nDim; ++m)
{
size_t idx = i * nDim + m;
// gradient of normalization
As.col(idx) -= ke.sign * nrm * nrm * dSurf[idx] * w * sfp * dPhi.dot(dPhi) * detJ; // d1/S * (psi + grad_psi) * grad_phi^2 * detJ
// gradient of stabilization term (surface)
As.col(idx) += ke.sign * nrm * w * dGradSfS[idx] * dPhi.dot(dPhi) * detJ; // 1/S * dgrad_psi * grad_phi^2 * detJ
// gradient of jacobian determinant
As.col(idx) += ke.sign * w * sfp * dPhi.dot(dPhi) * dDetJ[idx]; // psi * grad_phi^2 * ddetJ
As.col(idx) += ke.sign * nrm * w * sfp * dPhi.dot(dPhi) * dDetJ[idx]; // 1/S * (psi + grad_psi) * grad_phi^2 * ddetJ
}
}
}
......
dart/src/wKuttaResidual.h
View file @ 5ba134b1
......@@ -41,4 +41,4 @@ public:
};
} // namespace dart
#endif //WKUTTARESIDUAL_H
#endif // WKUTTARESIDUAL_H
dart/src/wLoadFunctional.h
View file @ 5ba134b1
......@@ -39,4 +39,4 @@ public:
};
} // namespace dart
#endif //WLOADFUNCTIONAL_H
#endif // WLOADFUNCTIONAL_H
dart/src/wNewton.cpp
View file @ 5ba134b1
......@@ -41,6 +41,7 @@
#include <tbb/global_control.h>
#include <tbb/parallel_for_each.h>
#include <tbb/spin_mutex.h>
#include <tbb/combinable.h>
#include <iomanip>
#include <deque>
......@@ -81,17 +82,18 @@ Newton::Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
*
* Solve the steady transonic Full Potential Equation
*/
bool Newton::run()
Status Newton::run()
{
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
// Display current freestream conditions
std::cout << std::setprecision(2)
<< "- Mach " << pbl->M_inf << ", "
<< pbl->alpha * 180 / 3.14159 << " deg AoA, "
<< pbl->beta * 180 / 3.14159 << " deg AoS"
<< std::endl;
if (verbose > 0)
std::cout << std::fixed << std::setprecision(2)
<< "Solving flow for Mach " << pbl->M_inf << ", "
<< pbl->alpha * 180 / 3.14159 << "deg AoA, "
<< pbl->beta * 180 / 3.14159 << "deg AoS"
<< std::endl;
// Initialize solver loop
nIt = 0; // iteration counter
......@@ -113,25 +115,28 @@ bool Newton::run()
// Linesearch
FpeLSFunction fpe(*this, phi, phi, dPhi_, rPhi_);
BankRoseLS ls(fpe);
//FleuryLS ls(fpe, 1., 1e-2, false);
// FleuryLS ls(fpe, 1., 1e-2, false);
ls.set(maxLsIt, lsTol, verbose);
// Display residual
std::cout << std::setw(8) << "N_Iter"
<< std::setw(8) << "L_Iter"
<< std::setw(8) << "f_eval"
<< std::setw(12) << "Cl"
<< std::setw(12) << "Cd"
<< std::setw(15) << "Rel_Res[phi]"
<< std::setw(15) << "Abs_Res[phi]" << std::endl;
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
<< std::setw(8) << 0
<< std::setw(8) << 1
<< std::setw(12) << "-"
<< std::setw(12) << "-"
<< std::setw(15) << log10(relRes)
<< std::setw(15) << log10(absRes) << "\n";
if (verbose > 0)
{
std::cout << std::setw(8) << "N_Iter"
<< std::setw(8) << "L_Iter"
<< std::setw(8) << "f_eval"
<< std::setw(12) << "Cl"
<< std::setw(12) << "Cd"
<< std::setw(15) << "Rel_Res[phi]"
<< std::setw(15) << "Abs_Res[phi]" << std::endl;
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
<< std::setw(8) << 0
<< std::setw(8) << 1
<< std::setw(12) << "-"
<< std::setw(12) << "-"
<< std::setw(15) << log10(relRes)
<< std::setw(15) << log10(absRes) << "\n";
}
do
{
......@@ -145,19 +150,19 @@ bool Newton::run()
{
case 0:
{
mCOv = 0.92;
mCOv = 0.925;
muCv = 2.;
break;
}
case 1:
{
mCOv = 0.935;
mCOv = 0.95;
muCv = 1.5;
break;
}
default:
{
mCOv = 0.95;
mCOv = 0.975;
muCv = 1.;
break;
}
......@@ -197,7 +202,7 @@ bool Newton::run()
Solver::computeLoad();
// Display residual (at each iteration)
if (verbose > 0)
if (verbose > 1)
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
......@@ -219,7 +224,7 @@ bool Newton::run()
} while (nIt < maxIt);
// Display residual (only last iteration)
if (verbose == 0)
if (verbose == 1)
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
......@@ -237,30 +242,36 @@ bool Newton::run()
// Check the solution
if (relRes < relTol || absRes < absTol)
{
std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 1)
if (verbose > 0)
std::cout << ANSI_COLOR_GREEN << "Newton solver converged." << ANSI_COLOR_RESET << std::endl;
if (verbose > 2)
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
return true;
if (verbose > 0)
std::cout << std::endl;
return Status::CONVERGED;
}
else if (std::isnan(relRes))
{
std::cout << ANSI_COLOR_RED << "Newton solver diverged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 1)
if (verbose > 0)
std::cout << ANSI_COLOR_RED << "Newton solver diverged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 2)
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
return false;
if (verbose > 0)
std::cout << std::endl;
return Status::FAILED;
}
else
{
std::cout << ANSI_COLOR_YELLOW << "Newton solver not fully converged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 1)
if (verbose > 0)
std::cout << ANSI_COLOR_YELLOW << "Newton solver not fully converged." << ANSI_COLOR_RESET << std::endl;
if (verbose > 2)
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
return true;
if (verbose > 0)
std::cout << std::endl;
return Status::MAXIT;
}
}
......@@ -271,11 +282,8 @@ bool Newton::run()
*/
void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
{
// Multithread
tbb::spin_mutex mutex;
// List of triplets to build Jacobian matrix
std::deque<Eigen::Triplet<double>> T;
// List of (thead-safe) triplets to build Jacobian matrix
tbb::combinable<std::deque<Eigen::Triplet<double>>> cmbT;
// Full Potential Equation with upwind bias: analytical derivatives
tms["0-Jbase"].start();
......@@ -283,25 +291,27 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
tbb::parallel_for_each(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
// Current element
Element *e = p.first;
//Upwind element
// Upwind element
Element *eU = p.second[0];
// Build elementary matrices
Eigen::MatrixXd Ae1, Ae2;
std::tie(Ae1, Ae2) = PotentialResidual::buildGradientFlow(*e, *eU, phi, *fluid, muCv, mCOv);
// Assembly (subsonic)
tbb::spin_mutex::scoped_lock lock(mutex);
std::deque<Eigen::Triplet<double>> &locT = cmbT.local();
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae1(ii, jj)));
locT.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae1(ii, jj)));
}
}
// Assembly (supersonic)
double mach = fluid->mach->eval(*e, phi, 0);
double machC = fluid->mach->eval(*e, phi, 0);
double machU = fluid->mach->eval(*eU, phi, 0);
double mach = (machC < machU) ? machU : machC;
if (mach > mCOv)
{
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
......@@ -310,7 +320,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
for (size_t jj = 0; jj < eU->nodes.size(); ++jj)
{
Node *nodj = eU->nodes[jj];
T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae2(ii, jj)));
locT.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae2(ii, jj)));
}
}
}
......@@ -328,19 +338,19 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
Eigen::MatrixXd Aue, Ale;
std::tie(Aue, Ale) = WakeResidual::buildGradientFlow(*we, phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
std::deque<Eigen::Triplet<double>> &locT = cmbT.local();
for (size_t ii = 0; ii < we->nRow; ++ii)
{
Node *nodi = we->wkN[ii].second;
for (size_t jj = 0; jj < we->nColUp; ++jj)
{
Node *nodj = we->volUpE->nodes[jj];
T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Aue(ii, jj)));
locT.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Aue(ii, jj)));
}
for (size_t jj = 0; jj < we->nColLw; ++jj)
{
Node *nodj = we->volLwE->nodes[jj];
T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, -Ale(ii, jj)));
locT.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, -Ale(ii, jj)));
}
}
});
......@@ -354,14 +364,14 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
// Build K
Eigen::MatrixXd Ae = KuttaResidual::buildGradientFlow(*ke, phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
std::deque<Eigen::Triplet<double>> &locT = cmbT.local();
for (size_t ii = 0; ii < ke->nRow; ++ii)
{
Node *nodi = ke->teN[ii].second;
for (size_t jj = 0; jj < ke->nCol; ++jj)
{
Node *nodj = ke->volE->nodes[jj];
T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae(ii, jj)));
locT.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae(ii, jj)));
}
}
});
......@@ -369,6 +379,10 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
tms["1-Jkutta"].stop();
// Build Jacobian matrix without BCs
std::deque<Eigen::Triplet<double>> T = cmbT.combine([&](std::deque<Eigen::Triplet<double>> x, std::deque<Eigen::Triplet<double>> const &y) {
x.insert(x.begin(), y.begin(), y.end());
return x;
});
J.setFromTriplets(T.begin(), T.end());
// Apply Dirichlet BCs
......@@ -392,7 +406,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
J.prune(0.);
J.makeCompressed();
if (verbose > 2)
if (verbose > 3)
std::cout << "J (" << J.rows() << "," << J.cols() << ") nnz=" << J.nonZeros() << "\n";
}
......@@ -502,7 +516,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
for (auto nod : dBC->nodes)
R(nod->row) = 0.;
if (verbose > 2)
if (verbose > 3)
std::cout << "R (" << R.size() << ")\n";
}
......
dart/src/wNewton.h
View file @ 5ba134b1
......@@ -49,7 +49,7 @@ private:
public:
Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol = 1e-6, double aTol = 1e-8, int mIt = 25, int nthrds = 1, int vrb = 1);
virtual bool run() override;
virtual Status run() override;
#ifndef SWIG
virtual void write(std::ostream &out) const override;
......@@ -63,4 +63,4 @@ private:
} // namespace dart
#endif //WNEWTON_H
#endif // WNEWTON_H
dart/src/wPicard.cpp
View file @ 5ba134b1
......@@ -75,17 +75,18 @@ Picard::Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
*
* Solve the steady subsonic Full Potential Equation
*/
bool Picard::run()
Status Picard::run()
{
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
// Display current freestream conditions
std::cout << std::setprecision(2)
<< "- Mach " << pbl->M_inf << ", "
<< pbl->alpha * 180 / 3.14159 << " deg AoA, "
<< pbl->beta * 180 / 3.14159 << " deg AoS"
<< std::endl;
if (verbose > 0)
std::cout << std::fixed << std::setprecision(2)
<< "Solving flow for Mach " << pbl->M_inf << ", "
<< pbl->alpha * 180 / 3.14159 << "deg AoA, "
<< pbl->beta * 180 / 3.14159 << "deg AoS"
<< std::endl;
// Initialize solver loop
nIt = 0;
......@@ -95,11 +96,12 @@ bool Picard::run()
Eigen::Map<Eigen::VectorXd> phi_(phi.data(), phi.size()), rPhi_(rPhi.data(), rPhi.size());
Eigen::VectorXd phiOld(phi.size());
std::cout << std::setw(8) << "N_Iter"
<< std::setw(12) << "Cl"
<< std::setw(12) << "Cd"
<< std::setw(15) << "Rel_Res[phi]"
<< std::setw(15) << "Abs_Res[phi]" << std::endl;
if (verbose > 0)
std::cout << std::setw(8) << "N_Iter"
<< std::setw(12) << "Cl"
<< std::setw(12) << "Cd"
<< std::setw(15) << "Rel_Res[phi]"
<< std::setw(15) << "Abs_Res[phi]" << std::endl;
do
{
......@@ -122,7 +124,7 @@ bool Picard::run()
Solver::computeLoad();
// Display residual (at each iteration)
if (verbose > 0 || nIt == 0)
if (verbose > 1 || (verbose == 1 && nIt == 0))
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
......@@ -152,7 +154,7 @@ bool Picard::run()
} while (nIt < maxIt);
// Display residual (only last iteration)
if (verbose == 0)
if (verbose == 1)
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
......@@ -168,36 +170,42 @@ bool Picard::run()
// Check the solution
if (relRes < relTol || absRes < absTol)
{
std::cout << ANSI_COLOR_GREEN << "Picard solver converged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 1)
if (verbose > 0)
std::cout << ANSI_COLOR_GREEN << "Picard solver converged." << ANSI_COLOR_RESET << std::endl;
if (verbose > 2)
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
return true;
if (verbose > 0)
std::cout << std::endl;
return Status::CONVERGED;
}
else if (std::isnan(relRes))
{
std::cout << ANSI_COLOR_RED << "Picard sovler diverged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 1)
if (verbose > 0)
std::cout << ANSI_COLOR_RED << "Picard sovler diverged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 2)
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
return false;
if (verbose > 0)
std::cout << std::endl;
return Status::FAILED;
}
else
{
std::cout << ANSI_COLOR_YELLOW << "Picard solver not fully converged!" << ANSI_COLOR_RESET << std::endl;
if (verbose > 1)
if (verbose > 0)
std::cout << ANSI_COLOR_YELLOW << "Picard solver not fully converged." << ANSI_COLOR_RESET << std::endl;
if (verbose > 2)
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
return true;
if (verbose > 0)
std::cout << std::endl;
return Status::MAXIT;
}
}
/**
* @brief Build LHS matrix and RHS vector for Picard iteration
*
*
* A = \int rho * grad_phi * grad_psi dV
* b = \int rho * grad_phi * n * psi dS
*/
......@@ -345,7 +353,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
A.prune(0.);
A.makeCompressed();
if (verbose > 2)
if (verbose > 3)
{
std::cout << "A (" << A.rows() << "," << A.cols() << ") nnz=" << A.nonZeros() << "\n";
std::cout << "b (" << b.size() << ")\n";
......
dart/src/wPicard.h
View file @ 5ba134b1
......@@ -40,7 +40,7 @@ public:
public:
Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol = 1e-6, double aTol = 1e-8, int mIt = 25, int nthrds = 1, int vrb = 1);
virtual bool run() override;
virtual Status run() override;
#ifndef SWIG
virtual void write(std::ostream &out) const override;
......@@ -52,4 +52,4 @@ private:
} // namespace dart
#endif //WPICARD_H
#endif // WPICARD_H
dart/src/wPotentialResidual.cpp
View file @ 5ba134b1
......@@ -64,7 +64,9 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
b += fluid.rho->eval(e, phi, k) * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
// Supersonic contribution
double mach = fluid.mach->eval(e, phi, 0);
double machC = fluid.mach->eval(e, phi, 0);
double machU = fluid.mach->eval(eU, phi, 0);
double mach = (machC < machU) ? machU : machC;
if (mach > mCO)
{
double mu = muC * (1 - (mCO * mCO) / (mach * mach)); // switching function
......@@ -110,7 +112,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
// Supersonic contribution
Eigen::MatrixXd A2;
double mach = fluid.mach->eval(e, phi, 0);
double machC = fluid.mach->eval(e, phi, 0);
double machU = fluid.mach->eval(eU, phi, 0);
double mach = (machC < machU) ? machU : machC;
if (mach > mCO)
{
// switching function and gradient
......@@ -133,11 +137,15 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
Eigen::VectorXd dPhi = e.computeGradient(phi, k);
// mu*drhoU*grad_phi*grad_psi
// mu * drhoU * grad_phi * grad_psi
A2 += mu * dRhoU * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj;
// mu*rhoU*grad_phi*grad_psi + dmu*(rhoU-rho)*grad_phi*grad_psi
// mu * rhoU * grad_phi * grad_psi
A1 += mu * rhoU * dSf.transpose() * dSf * wdj;
A1 += dmu * (rhoU - fluid.rho->eval(e, phi, k)) * fluid.mach->evalGrad(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj;
// dmu * (rhoU-rho) * grad_phi * grad_psi OR dmuU * (rhoU-rho) * grad_phi * grad_psi
if (machC >= machU)
A1 += dmu * (rhoU - fluid.rho->eval(e, phi, k)) * fluid.mach->evalGrad(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj;
else
A2 += dmu * (rhoU - fluid.rho->eval(e, phi, k)) * fluid.mach->evalGrad(eU, phi, k) * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj;
}
}
return std::make_tuple(A1, A2);
......@@ -181,7 +189,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
for (int m = 0; m < nDim; ++m)
{
size_t idx = i * nDim + m;
A1.col(idx) += w * dRho * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ; // drho * grad_phi*grad_psi * detj
A1.col(idx) += w * dRho * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ; // drho * grad_phi * grad_psi * detj
A1.col(idx) += w * rho * (iJ * dSf).transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // rho * dgrad_phi * grad_psi * detj
A1.col(idx) += w * rho * (-iJ * dJ[idx] * iJ * dSf).transpose() * dPhi * detJ; // rho * grad_phi * dgrad_psi * detj
A1.col(idx) += w * rho * (iJ * dSf).transpose() * dPhi * dDetJ[idx]; // rho * grad_phi * grad_psi * ddetj
......@@ -190,7 +198,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
}
// Supersonic contribution
Eigen::MatrixXd A2;
double mach = fluid.mach->eval(e, phi, 0);
double machC = fluid.mach->eval(e, phi, 0);
double machU = fluid.mach->eval(eU, phi, 0);
double mach = (machC < machU) ? machU : machC;
if (mach > mCO)
{
// switching function and gradient
......@@ -211,7 +221,6 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
// Density, Mach and Gradients
Eigen::VectorXd dPhi = e.computeGradient(phi, k);
double rho = fluid.rho->eval(e, phi, k);
double dM = fluid.mach->evalGrad(e, phi, k);
Eigen::MatrixXd const &dSf = cache.getDsf(k);
Eigen::MatrixXd const &iJ = e.getJinv(k);
// Jacobian gradients
......@@ -225,10 +234,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
for (size_t i = 0; i < e.nodes.size(); ++i)
{
size_t idx = i * nDim + m;
A1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // mu * rhoU * dgrad_phi * grad_psi * detj
A1.col(idx) += w * mu * rhoU * (-iJ * dJ[idx] * iJ * dSf).transpose() * dPhi * detJ; // mu * rhoU * grad_phi * dgrad_psi * detj
A1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * dPhi * dDetJ[idx]; // mu * rhoU * grad_phi * grad_psi * ddetj
A1.col(idx) += w * (-rho + rhoU) * dmu * dM * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ; // (rhoU-rho) * dmu * grad_phi * grad_psi * detj
A1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // mu * rhoU * dgrad_phi * grad_psi * detj
A1.col(idx) += w * mu * rhoU * (-iJ * dJ[idx] * iJ * dSf).transpose() * dPhi * detJ; // mu * rhoU * grad_phi * dgrad_psi * detj
A1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * dPhi * dDetJ[idx]; // mu * rhoU * grad_phi * grad_psi * ddetj
}
for (size_t i = 0; i < eU.nodes.size(); ++i)
{
......@@ -236,6 +244,21 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
A2.col(idx) += w * mu * dRhoU * dPhiU.transpose() * (-iJU * dJU[idx] * dPhiU) * (iJ * dSf).transpose() * dPhi * detJ; // mu * drhoU * grad_phi * grad_psi * detj
}
}
// Gradient of swithcing function, (rhoU-rho) * dmu * grad_phi * grad_psi * detj OR (rhoU-rho) * dmuU * grad_phi * grad_psi * detj
if (machC >= machU)
{
double dM = fluid.mach->evalGrad(e, phi, k);
for (int m = 0; m < nDim; ++m)
for (size_t i = 0; i < e.nodes.size(); ++i)
A1.col(i * nDim + m) += w * (-rho + rhoU) * dmu * dM * dPhi.transpose() * (-iJ * dJ[i * nDim + m] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ;
}
else
{
double dM = fluid.mach->evalGrad(eU, phi, k);
for (int m = 0; m < nDim; ++m)
for (size_t i = 0; i < eU.nodes.size(); ++i)
A2.col(i * nDim + m) += w * (-rho + rhoU) * dmu * dM * dPhiU.transpose() * (-iJU * dJU[i * nDim + m] * dPhiU) * (iJ * dSf).transpose() * dPhi * detJ;
}
}
}
return std::make_tuple(A1, A2);
......
dart/src/wPotentialResidual.h
View file @ 5ba134b1
......@@ -43,4 +43,4 @@ public:
};
} // namespace dart
#endif //WPOTENTIALRESIDUAL_H
#endif // WPOTENTIALRESIDUAL_H
dart/src/wProblem.cpp
View file @ 5ba134b1
......@@ -194,63 +194,42 @@ void Problem::check() const
// Sanity checks
if (fluid == NULL)
throw std::runtime_error("No Fluid provided!\n");
if (bodies.empty())
throw std::runtime_error("No Body provided!\n");
if (iIC == NULL)
throw std::runtime_error("No Initial Condition provided!\n");
if (fBCs.empty())
throw std::runtime_error("No Freestream B.C. provided!\n");
if (wBCs.empty())
throw std::runtime_error("No Wake B.C. provided!\n");
if (kSCs.empty())
throw std::runtime_error("No Kutta condition provided!\n");
// Three-dimension problem
if (nDim == 3)
{
// Volume element type
for (auto e : fluid->tag->elems)
{
if (e->type() != ELTYPE::TETRA4)
{
std::stringstream err;
err << "3D solver is only implemented for volume elements of type "
<< ELTYPE::TETRA4 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
if (e->type() != ElType::TETRA4)
throwUnsupElType(e, "3", "volume", ElType::TETRA4);
// Surface element type (Body)
for (auto surf : bodies)
for (auto e : surf->groups[0]->tag->elems)
if (e->type() != ElType::TRI3)
throwUnsupElType(e, "3", "surface", ElType::TRI3);
// Surface element type (Freestream B.C.)
for (auto surf : fBCs)
{
for (auto e : surf->tag->elems)
{
if (e->type() != ELTYPE::TRI3)
{
std::stringstream err;
err << "3D solver is only implemented for surface elements of type "
<< ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
}
if (e->type() != ElType::TRI3)
throwUnsupElType(e, "3", "surface", ElType::TRI3);
// Surface element type (Wake B.C.)
for (auto surf : wBCs)
{
for (auto e : surf->groups[0]->tag->elems)
{
if (e->type() != ELTYPE::TRI3)
{
std::stringstream err;
err << "3D solver is only implemented for surface elements of type "
<< ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
if (e->type() != ElType::TRI3)
throwUnsupElType(e, "3", "surface", ElType::TRI3);
for (auto e : surf->groups[1]->tag->elems)
{
if (e->type() != ELTYPE::TRI3)
{
std::stringstream err;
err << "3D solver is only implemented for surface elements of type "
<< ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
if (e->type() != ElType::TRI3)
throwUnsupElType(e, "3", "surface", ElType::TRI3);
}
// Blowing B.C.
if (!bBCs.empty())
......@@ -261,52 +240,27 @@ void Problem::check() const
{
// Volume element type
for (auto e : fluid->tag->elems)
{
if (e->type() != ELTYPE::TRI3)
{
std::stringstream err;
err << "2D solver is only implemented for volume elements of type "
<< ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
if (e->type() != ElType::TRI3)
throwUnsupElType(e, "2", "volume", ElType::TRI3);
// Surface element type (Body)
for (auto surf : bodies)
for (auto e : surf->groups[0]->tag->elems)
if (e->type() != ElType::LINE2)
throwUnsupElType(e, "2", "surface", ElType::LINE2);
// Surface element type (Freestream B.C.)
for (auto surf : fBCs)
{
for (auto e : surf->tag->elems)
{
if (e->type() != ELTYPE::LINE2)
{
std::stringstream err;
err << "2D solver is only implemented for surface elements of type "
<< ELTYPE::LINE2 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
}
if (e->type() != ElType::LINE2)
throwUnsupElType(e, "2", "surface", ElType::LINE2);
// Surface element type (Wake B.C.)
for (auto surf : wBCs)
{
for (auto e : surf->groups[0]->tag->elems)
{
if (e->type() != ELTYPE::LINE2)
{
std::stringstream err;
err << "2D solver is only implemented for surface elements of type "
<< ELTYPE::LINE2 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
if (e->type() != ElType::LINE2)
throwUnsupElType(e, "2", "surface", ElType::LINE2);
for (auto e : surf->groups[1]->tag->elems)
{
if (e->type() != ELTYPE::LINE2)
{
std::stringstream err;
err << "2D solver is only implemented for surface elements of type "
<< ELTYPE::LINE2 << " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
}
if (e->type() != ElType::LINE2)
throwUnsupElType(e, "2", "surface", ElType::LINE2);
}
}
else
......@@ -318,6 +272,15 @@ void Problem::check() const
}
}
void Problem::throwUnsupElType(Element const *e, std::string const &pdim, std::string const &edim, ElType type) const
{
std::stringstream err;
err << pdim << "D solver is only implemented for "
<< edim << " elements of type " << type
<< " (" << e->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
void Problem::write(std::ostream &out) const
{
out << "dart::Problem parameters"
......
dart/src/wProblem.h
View file @ 5ba134b1
......@@ -19,6 +19,7 @@
#include "dart.h"
#include "wObject.h"
#include "wElement.h"
#include "wF1Ct.h"
#include <memory>
#include <iostream>
......@@ -87,8 +88,11 @@ public:
void initGradElems();
virtual void write(std::ostream &out) const override;
#endif
private:
void throwUnsupElType(tbox::Element const *e, std::string const &pdim, std::string const &edim, tbox::ElType type) const;
};
} // namespace dart
#endif //WPROBLEM_H
#endif // WPROBLEM_H
dart/src/wSolver.cpp
View file @ 5ba134b1
......@@ -21,6 +21,7 @@
#include "wAssign.h"
#include "wBody.h"
#include "wWake.h"
#include "wKutta.h"
#include "wLoadFunctional.h"
#include "wMshData.h"
......@@ -35,6 +36,7 @@
#include <tbb/spin_mutex.h>
#include <iomanip>
#include <fstream>
using namespace tbox;
using namespace dart;
......@@ -63,9 +65,9 @@ Solver::Solver(std::shared_ptr<Problem> _pbl,
<< "** _ /_/ /_ ___ | __ _/_ / **\n"
<< "** /_____/ /_/ |_/_/ |_| /_/ **\n"
<< "*******************************************************************************\n"
<< "** Hi! My name is DART v1.1.0-21.09 **\n"
<< "** Hi! My name is DART v1.3.0-24.12 **\n"
<< "** Adrien Crovato **\n"
<< "** ULiege 2018-2021 **\n"
<< "** ULiege 2018-2024 **\n"
<< "*******************************************************************************\n"
<< std::endl;
......@@ -110,6 +112,8 @@ Solver::Solver(std::shared_ptr<Problem> _pbl,
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"
......@@ -150,7 +154,7 @@ void Solver::reset()
/**
* @brief Dummy run
*/
bool Solver::run()
Status Solver::run()
{
throw std::runtime_error("dart::Solver::run not implemented!\n");
}
......@@ -158,11 +162,16 @@ bool Solver::run()
/**
* @brief Write the results
*/
void Solver::save(MshExport *mshWriter, int n)
void Solver::save(MshExport *mshWriter, std::string const &suffix)
{
// Write files
std::cout << "Saving files... " << std::endl;
// setup results
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;
......@@ -171,20 +180,30 @@ void Solver::save(MshExport *mshWriter, int n)
results.scalars_at_nodes["rho"] = &rho;
results.scalars_at_nodes["Mach"] = &M;
results.scalars_at_nodes["Cp"] = &Cp;
// save (whole mesh and bodies)
if (n > 0)
{
mshWriter->save(pbl->msh->name + "_" + std::to_string(n), results);
for (auto bnd : pbl->bodies)
bnd->save(bnd->groups[0]->tag->name + "_" + std::to_string(n), results);
}
else
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)
{
mshWriter->save(pbl->msh->name, results);
for (auto bnd : pbl->bodies)
bnd->save(bnd->groups[0]->tag->name, results);
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";
}
std::cout << std::endl;
outfile.close();
std::cout << "done." << std::endl;
}
/**
......@@ -276,6 +295,6 @@ void Solver::computeFlow()
});
// 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;
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;
}
dart/src/wSolver.h
View file @ 5ba134b1
......@@ -29,6 +29,16 @@
namespace dart
{
/**
* @brief Exit code of DART solver
*/
enum class Status
{
CONVERGED = 0, // fully converged to prescribed (relative or absolute) tolerance
MAXIT = 1, // not fully converged (max. number of iterations exceeded)
FAILED = 2 // NaN in the solution vector
};
/**
* @brief Base class for full potential solvers
* @authors Adrien Crovato
......@@ -67,8 +77,10 @@ public:
virtual ~Solver();
void reset();
virtual bool run();
void save(tbox::MshExport *mshWriter, int n = 0);
virtual Status run();
void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
int getRow(size_t i) const { return rows[i]; };
protected:
void computeFlow();
......@@ -76,4 +88,4 @@ protected:
};
} // namespace dart
#endif //WSOLVER_H
#endif // WSOLVER_H
dart/src/wWake.cpp
View file @ 5ba134b1
......@@ -112,10 +112,10 @@ void Wake::connectNodes()
// Create the node map
for (auto nUp : nodesUp)
{
//std::cout << "MASTER: " << nUp->pos(0) << ',' << nUp->pos(1) << ',' << nUp->pos(2) << std::endl;
// std::cout << "MASTER: " << nUp->pos(0) << ',' << nUp->pos(1) << ',' << nUp->pos(2) << std::endl;
for (auto nLw : nodesLw)
{
//std::cout << " SLAVE: " << nLw->pos(0) << ',' << nLw->pos(1) << ',' << nLw->pos(2) << std::endl;
// std::cout << " SLAVE: " << nLw->pos(0) << ',' << nLw->pos(1) << ',' << nLw->pos(2) << std::endl;
if ((nUp->pos - nLw->pos).norm() <= 1e-14)
{
nodMap[nUp] = nLw;
......
dart/src/wWake.h
View file @ 5ba134b1
......@@ -59,4 +59,4 @@ private:
} // namespace dart
#endif //WWAKE_H
#endif // WWAKE_H
dart/src/wWakeElement.cpp
View file @ 5ba134b1
......@@ -28,35 +28,35 @@ WakeElement::WakeElement(size_t _no, Element *&_surUpE, Element *&_surLwE,
// Sanity checks
if (surUpE->type() != surLwE->type())
throw std::runtime_error("WakeElement: upper and lower surface element are not the same type!\n");
if (surUpE->type() == ELTYPE::LINE2)
if (surUpE->type() == ElType::LINE2)
nDim = 2;
else if (surUpE->type() == ELTYPE::TRI3)
else if (surUpE->type() == ElType::TRI3)
nDim = 3;
else
{
std::stringstream err;
err << "WakeElement only implemented for surface elements of type "
<< ELTYPE::LINE2 << " or " << ELTYPE::TRI3
<< ElType::LINE2 << " or " << ElType::TRI3
<< " (" << surUpE->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
if (volUpE->type() == ELTYPE::TRI3 || volUpE->type() == ELTYPE::TETRA4)
if (volUpE->type() == ElType::TRI3 || volUpE->type() == ElType::TETRA4)
nColUp = volUpE->nodes.size();
else
{
std::stringstream err;
err << "WakeElement only implemented for volume elements of type "
<< ELTYPE::TRI3 << " or " << ELTYPE::TETRA4
<< ElType::TRI3 << " or " << ElType::TETRA4
<< " (" << volUpE->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
if (volLwE->type() == ELTYPE::TRI3 || volLwE->type() == ELTYPE::TETRA4)
if (volLwE->type() == ElType::TRI3 || volLwE->type() == ElType::TETRA4)
nColLw = volLwE->nodes.size();
else
{
std::stringstream err;
err << "WakeElement only implemented for volume elements of type "
<< ELTYPE::TRI3 << " or " << ELTYPE::TETRA4
<< ElType::TRI3 << " or " << ElType::TETRA4
<< " (" << volLwE->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
......
dart/src/wWakeElement.h
View file @ 5ba134b1
......@@ -30,8 +30,8 @@ namespace dart
* @brief Super-element to apply wake boundary conditions
*
* Made up of two surface elements on the wake and the connected volume
* elements. Only Line2 (surface) and Tri3 (volume) elements are curently
* implemented.
* elements. Only Line2/Tri3 (surface) and Tri3/Tetra4 (volume) elements are
* currently implemented.
* Since Tri3 or Tetra4 are used, the evalutaion of the volume integrals is
* performed at dummy Gauss point #0. If volume elements with non-constant shape
* function derivatives are used, those derivates should be recomputed at the
......@@ -67,4 +67,4 @@ private:
} // namespace dart
#endif //WWAKEELEMENT_H
#endif // WWAKEELEMENT_H
dart/src/wWakeResidual.cpp
View file @ 5ba134b1
......@@ -66,7 +66,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
sfp(i) = dSfp(i) + sf(we.wkN[i].first);
sfp(i) = sf(we.wkN[i].first) + dSfp(i);
// wake contribution
Au += sfp * Aup * sGauss.getW(k) * we.surUpE->getDetJ(k);
Al += sfp * Alw * sGauss.getW(k) * we.surUpE->getDetJ(k);
......@@ -99,10 +99,8 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
Eigen::VectorXd dSfp(nRow);
dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
// Velocity squared
Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0);
Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0);
double dPhi2u = dPhiU.squaredNorm();
double dPhi2l = dPhiL.squaredNorm();
double dPhi2u = we.volUpE->computeGradient(phi, 0).squaredNorm();
double dPhi2l = we.volLwE->computeGradient(phi, 0).squaredNorm();
// Build
Eigen::VectorXd bu = Eigen::VectorXd::Zero(nRow);
......@@ -113,7 +111,7 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
sfp(i) = dSfp(i) + sf(we.wkN[i].first);
sfp(i) = sf(we.wkN[i].first) + dSfp(i);
// wake contribution
bu += sfp * dPhi2u * sGauss.getW(k) * we.surUpE->getDetJ(k);
bl += sfp * dPhi2l * sGauss.getW(k) * we.surUpE->getDetJ(k);
......