Skip to content
Snippets Groups Projects
Commit 3e08876f authored by Adrien Crovato's avatar Adrien Crovato
Browse files

Remove Kutta stabilization for 2D

parent 4d61def3
No related branches found
No related tags found
No related merge requests found
...@@ -42,13 +42,16 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do ...@@ -42,13 +42,16 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
size_t nRow = ke.nRow; size_t nRow = ke.nRow;
size_t nDim = ke.nDim; size_t nDim = ke.nDim;
// Stabilization // Stabilization (disabled in 2D)
Eigen::MatrixXd dSf(nRow, nDim); Eigen::VectorXd dSfp = Eigen::VectorXd::Zero(nRow);
for (size_t i = 0; i < nRow; ++i) if (nDim == 3)
for (size_t j = 0; j < nDim; ++j) {
dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first)); Eigen::MatrixXd dSf(nRow, nDim);
Eigen::VectorXd dSfp(nRow); for (size_t i = 0; i < nRow; ++i)
dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim); 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 // Velocity
Eigen::RowVectorXd Aup = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).transpose() * iJ * dSfV; Eigen::RowVectorXd Aup = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).transpose() * iJ * dSfV;
...@@ -85,13 +88,16 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double> ...@@ -85,13 +88,16 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
size_t nRow = ke.nRow; size_t nRow = ke.nRow;
size_t nDim = ke.nDim; size_t nDim = ke.nDim;
// Stabilization // Stabilization (disabled in 2D)
Eigen::MatrixXd dSf(nRow, nDim); Eigen::VectorXd dSfp = Eigen::VectorXd::Zero(nRow);
for (size_t i = 0; i < nRow; ++i) if (nDim == 3)
for (size_t j = 0; j < nDim; ++j) {
dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first)); Eigen::MatrixXd dSf(nRow, nDim);
Eigen::VectorXd dSfp(nRow); for (size_t i = 0; i < nRow; ++i)
dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim); 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 // Velocity
double dPhi2 = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).squaredNorm(); double dPhi2 = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).squaredNorm();
...@@ -151,31 +157,37 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku ...@@ -151,31 +157,37 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
size_t nCol = ke.nCol; size_t nCol = ke.nCol;
size_t nSur = ke.surE->nodes.size(); size_t nSur = ke.surE->nodes.size();
// Stabilization term // Stabilization (disabled in 2D)
Eigen::VectorXd vi = Eigen::Vector3d(1, 0, 0).middleRows(0, nDim); Eigen::VectorXd gradSfp = Eigen::VectorXd::Zero(nRow);
double sqSurf = sqrt(surf); std::vector<Eigen::VectorXd> dGradSfV(nDim * nCol, Eigen::VectorXd::Zero(nRow));
Eigen::MatrixXd dSf(nRow, nDim); std::vector<Eigen::VectorXd> dGradSfS(nDim * nSur, Eigen::VectorXd::Zero(nRow));
for (size_t i = 0; i < nRow; ++i) if (nDim == 3)
for (size_t m = 0; m < nDim; ++m)
dSf(i, m) = dSfV(m, ke.vsMap.at(ke.teN[i].first));
Eigen::VectorXd gradSfp = 0.5 * sqSurf * dSf * iJ.transpose() * vi;
// Gradient of stabilization term (volume and surface)
std::vector<Eigen::VectorXd> dGradSfV(nDim * nCol);
for (size_t i = 0; i < nCol; ++i)
{ {
for (size_t m = 0; m < nDim; ++m) // 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)
{ {
size_t idx = i * nDim + m; for (size_t m = 0; m < nDim; ++m)
dGradSfV[idx] = 0.5 * sqSurf * dSf * (-iJ * dJ[idx] * iJ).transpose() * vi; // S * dgrad_psi {
size_t idx = i * nDim + m;
dGradSfV[idx] = 0.5 * sqSurf * dSf * (-iJ * dJ[idx] * iJ).transpose() * vi; // S * dgrad_psi
}
} }
} // surface gradient
std::vector<Eigen::VectorXd> dGradSfS(nDim * nSur); for (size_t i = 0; i < nSur; ++i)
for (size_t i = 0; i < nSur; ++i)
{
for (size_t m = 0; m < nDim; ++m)
{ {
size_t idx = i * nDim + m; for (size_t m = 0; m < nDim; ++m)
dGradSfS[idx] = 0.5 * 0.5 / sqSurf * dSurf[idx] * dSf * iJ.transpose() * vi; // dS * grad_psi {
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 // Velocity and normalization
......
...@@ -114,17 +114,17 @@ def main(): ...@@ -114,17 +114,17 @@ def main():
if M_inf == 0. and alpha == 0*np.pi/180: if M_inf == 0. and alpha == 0*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 6.9, 1e-2)) tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 6.9, 1e-2))
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.0, 1e-3)) tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.0, 1e-3))
tests.add(CTest('dCl_dMsh', dClX, 55.1, 1e-2)) tests.add(CTest('dCl_dMsh', dClX, 36.9, 1e-2))
tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-2)) tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-2))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, 7.1, 1e-2)) # FD 7.1364 (step = 1e-6) tests.add(CTest('dCl_dAoA (msh)', dClAoA, 6.9, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.0, 1e-3)) tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.0, 1e-3))
elif M_inf == 0.7 and alpha == 2*np.pi/180: elif M_inf == 0.7 and alpha == 2*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 11.8, 1e-2)) # FD 11.835 (step = 1e-5) tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 11.8, 1e-2)) # FD 11.835 (step = 1e-5)
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.127, 1e-2)) # 0.12692 (step = 1e-5) tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.127, 1e-2)) # 0.12692 (step = 1e-5)
tests.add(CTest('dCl_dMsh', dClX, 83.6, 1e-2)) tests.add(CTest('dCl_dMsh', dClX, 61.2, 1e-2))
tests.add(CTest('dCd_dMsh', dCdX, 0.904, 1e-2)) tests.add(CTest('dCd_dMsh', dCdX, 0.830, 1e-2))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, 12.1, 1e-2)) # FD 12.1443 (step = 1e-6) tests.add(CTest('dCl_dAoA (msh)', dClAoA, 11.8, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.129, 1e-2)) # FD 0.1289 (step = 1e-6) tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.127, 1e-2))
else: else:
raise Exception('Test not defined for this flow') raise Exception('Test not defined for this flow')
tests.run() tests.run()
......
...@@ -39,7 +39,6 @@ import dart.utils as floU ...@@ -39,7 +39,6 @@ import dart.utils as floU
import dart.default as floD import dart.default as floD
import dart.viscous.solver as floVS import dart.viscous.solver as floVS
import dart.viscous.coupler as floVC import dart.viscous.coupler as floVC
import tbox
import tbox.utils as tboxU import tbox.utils as tboxU
import fwk import fwk
from fwk.testing import * from fwk.testing import *
...@@ -53,11 +52,10 @@ def main(): ...@@ -53,11 +52,10 @@ def main():
# define flow variables # define flow variables
Re = 1e7 Re = 1e7
alpha = 5*math.pi/180 alpha = 5*math.pi/180
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
M_inf = 0.5 M_inf = 0.5
c_ref = 1 c_ref = 1
dim = 2 dim = 2
tol = 1e-3 # tolerance for the VII (usually one drag count) #TODO tolerance has been decreased with new Kutta, hopefully will increase with Paul's dev tol = 2e-3 # tolerance for the VII (usually one drag count) #TODO tolerance has been decreased with new Kutta, hopefully will increase with Paul's dev
# mesh the geometry # mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
......
...@@ -44,7 +44,7 @@ class Coupler: ...@@ -44,7 +44,7 @@ class Coupler:
# run inviscid solver # run inviscid solver
self.isolver.run() self.isolver.run()
print('--- Viscous solver parameters ---') print('--- Viscous solver parameters ---')
print('Tolerance (drag count):', self.tol * 1000) print('Tolerance (drag count):', self.tol * 1e4)
print('--- Viscous problem definition ---') print('--- Viscous problem definition ---')
print('Reynolds number:', self.vsolver.Re) print('Reynolds number:', self.vsolver.Re)
print('') print('')
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment