diff --git a/dart/src/wKuttaResidual.cpp b/dart/src/wKuttaResidual.cpp
index 6d0b9e3b806bd3ca058e16499b9adc56284465db..8884b5f7b7aa76db48ada9e9107c9df13935f608 100644
--- a/dart/src/wKuttaResidual.cpp
+++ b/dart/src/wKuttaResidual.cpp
@@ -42,13 +42,16 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
size_t nRow = ke.nRow;
size_t nDim = ke.nDim;
- // Stabilization
- 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));
- Eigen::VectorXd dSfp(nRow);
- dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, 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 = 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>
size_t nRow = ke.nRow;
size_t nDim = ke.nDim;
- // Stabilization
- 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));
- Eigen::VectorXd dSfp(nRow);
- dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, 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
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
size_t nCol = ke.nCol;
size_t nSur = ke.surE->nodes.size();
- // Stabilization 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));
- 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)
+ // 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)
{
- 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;
- dGradSfV[idx] = 0.5 * sqSurf * dSf * (-iJ * dJ[idx] * iJ).transpose() * vi; // S * dgrad_psi
+ 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
+ }
}
- }
- std::vector<Eigen::VectorXd> dGradSfS(nDim * nSur);
- for (size_t i = 0; i < nSur; ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
+ // surface gradient
+ for (size_t i = 0; i < nSur; ++i)
{
- size_t idx = i * nDim + m;
- dGradSfS[idx] = 0.5 * 0.5 / sqSurf * dSurf[idx] * dSf * iJ.transpose() * vi; // dS * grad_psi
+ 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
diff --git a/dart/tests/adjoint.py b/dart/tests/adjoint.py
index 6d35b848abe2b85bef09d01625878dbfce1115f8..5788fe27441057128a222cd330e4d421add55673 100644
--- a/dart/tests/adjoint.py
+++ b/dart/tests/adjoint.py
@@ -114,17 +114,17 @@ def main():
if M_inf == 0. and alpha == 0*np.pi/180:
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('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('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))
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('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('dCd_dMsh', dCdX, 0.904, 1e-2))
- tests.add(CTest('dCl_dAoA (msh)', dClAoA, 12.1, 1e-2)) # FD 12.1443 (step = 1e-6)
- tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.129, 1e-2)) # FD 0.1289 (step = 1e-6)
+ tests.add(CTest('dCl_dMsh', dClX, 61.2, 1e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.830, 1e-2))
+ tests.add(CTest('dCl_dAoA (msh)', dClAoA, 11.8, 1e-2))
+ tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.127, 1e-2))
else:
raise Exception('Test not defined for this flow')
tests.run()
diff --git a/dart/tests/bli.py b/dart/tests/bli.py
index 9014bcc0ab9feac43ed1dcb451d597950fd2cf60..ae4d0c9ec5ce1d7ce9888b47678003cde07d046c 100644
--- a/dart/tests/bli.py
+++ b/dart/tests/bli.py
@@ -39,7 +39,6 @@ import dart.utils as floU
import dart.default as floD
import dart.viscous.solver as floVS
import dart.viscous.coupler as floVC
-import tbox
import tbox.utils as tboxU
import fwk
from fwk.testing import *
@@ -53,11 +52,10 @@ def main():
# define flow variables
Re = 1e7
alpha = 5*math.pi/180
- U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
M_inf = 0.5
c_ref = 1
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
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
diff --git a/dart/viscous/coupler.py b/dart/viscous/coupler.py
index 269348dad08f5fb9362d2f2e05a7ce0c355e5585..6c339cfa5595818379c30cd17cfda7b66b593e4d 100644
--- a/dart/viscous/coupler.py
+++ b/dart/viscous/coupler.py
@@ -44,7 +44,7 @@ class Coupler:
# run inviscid solver
self.isolver.run()
print('--- Viscous solver parameters ---')
- print('Tolerance (drag count):', self.tol * 1000)
+ print('Tolerance (drag count):', self.tol * 1e4)
print('--- Viscous problem definition ---')
print('Reynolds number:', self.vsolver.Re)
print('')