diff --git a/sdpm/api/core.py b/sdpm/api/core.py
index d7d6b38764efba6dad38097667848bda0192dd1b..f0bbb3c07a8e7cea64740e3e0db8b3ae10a12b71 100644
--- a/sdpm/api/core.py
+++ b/sdpm/api/core.py
@@ -50,7 +50,7 @@ def init_sdpm(cfg, use_ad=False):
# Imports
import sdpm
from sdpm.gmsh import MeshLoader
- from sdpm.utils import interpolate_modes
+ from sdpm.utils import interpolate_modes, interpolate_pressure
import math
# Sanity checks
@@ -79,11 +79,14 @@ def init_sdpm(cfg, use_ad=False):
x_modes, amp_modes = interpolate_modes(cfg['Modes'], _bdy.getNodes())
for i in range(cfg['Num_modes']):
_bdy.setFlexibleMotion(i, 1 / amp_modes[i], x_modes[:, 3*i:3*i+3])
+ if 'Transonic_pressure_grad' in cfg:
+ dcp = interpolate_pressure(cfg['Transonic_pressure_grad'], _bdy.getNodes())
+ _bdy.setTransonicPressureGrad(dcp[:, 0])
_pbl.addBody(_bdy)
# Solver
vrb = cfg['Verb_lvl'] if 'Verb_lvl' in cfg else 1
- _sol = sdpm.Solver(_pbl, vrb)
+ _sol = sdpm.Solver(_pbl, vrb) if 'Transonic_pressure_grad' not in cfg else sdpm.SolverTransonic(_pbl, vrb)
_adj = sdpm.Adjoint(_sol) if use_ad else None
# Return
diff --git a/sdpm/src/sdpmSolverTransonic.cpp b/sdpm/src/sdpmSolverTransonic.cpp
index 7a1512249872cbee3ca207adfb0048c4ed9a1750..87835127bd38bfb02d827383a3f8ff974560ad69 100644
--- a/sdpm/src/sdpmSolverTransonic.cpp
+++ b/sdpm/src/sdpmSolverTransonic.cpp
@@ -22,7 +22,7 @@
using namespace sdpm;
-SolverTransonic::SolverTransonic(Problem &pbl) : Solver(pbl), _dCorr(sdpmVectorXd::Ones(_nP))
+SolverTransonic::SolverTransonic(Problem &pbl, int vrb) : Solver(pbl, vrb), _dCorr(sdpmVectorXd::Ones(_nP))
{
// Build node row id
int i = 0;
@@ -40,9 +40,10 @@ void SolverTransonic::run()
computeCorrection();
std::cout << "done." << std::endl;
- std::cout << "Running transonic SDPM solver...\n"
- << "-- Mean flow" << std::endl;
- std::cout << "computing sources... " << std::flush;
+ std::cout << "Running transonic SDPM solver..." << std::flush;
+ if (_vrb > 0)
+ std::cout << "\n-- Mean flow\n"
+ << "computing sources... " << std::flush;
sdpmDouble beta = _pbl.getCompressibility();
sdpmVector3d ui = _pbl.getDragDir();
sdpmVectorXd etauxy(_nP);
@@ -55,23 +56,27 @@ void SolverTransonic::run()
etauz(_rows[e]) = ui[2] * e->getCompressibleNormal()[2];
}
}
- std::cout << "done." << std::endl;
- std::cout << "solving for doublets... " << std::flush;
+ if (_vrb > 0)
+ std::cout << "done.\n"
+ << "solving for doublets... " << std::flush;
sdpmVectorXd rhsxy = _B0 * etauxy;
sdpmVectorXd rhsz = _B0 * etauz;
sdpmVectorXd emuxy = solve(_A0, rhsxy);
sdpmVectorXd emuz = solve(_A0, rhsz);
- std::cout << "done." << std::endl;
- std::cout << "computing flow and loads on bodies... " << std::flush;
+ if (_vrb > 0)
+ std::cout << "done.\n"
+ << "computing flow and loads on bodies... " << std::flush;
post(etauxy + etauz, emuxy + emuz.cwiseProduct(_dCorr));
- std::cout << "done." << std::endl;
+ if (_vrb > 0)
+ std::cout << "done." << std::endl;
for (size_t imd = 0; imd < _pbl.getNumModes(); ++imd)
{
for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
{
- std::cout << "-- Mode #" << imd << ", frequency #" << ifq << std::endl;
- std::cout << "computing sources... " << std::flush;
+ if (_vrb > 0)
+ std::cout << "-- Mode #" << imd << ", frequency #" << ifq << "\n"
+ << "computing sources... " << std::flush;
sdpmVectorXcd etau1(_nP);
sdpmVectorXcd etau1xy(_nP);
sdpmVectorXcd etau1z(_nP);
@@ -88,9 +93,10 @@ void SolverTransonic::run()
etau1z(_rows[elems[i]]) = (uc0[i](2) + uc1[i](2) * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]) * elems[i]->getCompressibleNormal()(2);
}
}
- std::cout << "done." << std::endl;
+ if (_vrb > 0)
+ std::cout << "done.\n"
+ << "solving for doublets... " << std::flush;
// Split complex set of equations into real set (twice the size)
- std::cout << "solving for doublets... " << std::flush;
sdpmMatrixXd A(2 * _nP, 2 * _nP);
A.block(0, 0, _nP, _nP) = _A[ifq].real();
A.block(_nP, _nP, _nP, _nP) = _A[ifq].real();
@@ -106,22 +112,29 @@ void SolverTransonic::run()
sdpmVectorXd muz = solve(A, rhs1z);
sdpmVectorXcd emu1xy = muxy.block(0, 0, _nP, 1) + sdpmComplex(0, 1) * muxy.block(_nP, 0, _nP, 1);
sdpmVectorXcd emu1z = muz.block(0, 0, _nP, 1) + sdpmComplex(0, 1) * muz.block(_nP, 0, _nP, 1);
- std::cout << "done." << std::endl;
- std::cout << "computing flow and loads on bodies... " << std::flush;
+ if (_vrb > 0)
+ std::cout << "done.\n"
+ << "computing flow and loads on bodies... " << std::flush;
sdpmVectorXcd mu1c = emu1xy + _dCorr.cwiseProduct(emu1z);
post(imd, ifq, etau1, mu1c);
- std::cout << "done." << std::endl;
+ if (_vrb > 0)
+ std::cout << "done." << std::endl;
}
}
// GAF
- std::cout << "-- Generalized Aerodynamic Force matrices" << std::endl;
+ if (_vrb > 0)
+ std::cout << "-- Generalized Aerodynamic Force matrices" << std::endl;
for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
{
- std::cout << "computing GAF at frequency #" << ifq << "... " << std::flush;
+ if (_vrb > 0)
+ std::cout << "computing GAF at frequency #" << ifq << "... " << std::flush;
computeGAF(ifq);
- std::cout << "done." << std::endl;
+ if (_vrb > 0)
+ std::cout << "done." << std::endl;
}
+ if (_vrb == 0)
+ std::cout << " done.";
std::cout << std::endl;
}
diff --git a/sdpm/src/sdpmSolverTransonic.h b/sdpm/src/sdpmSolverTransonic.h
index 7dcbea0b531ac8da3fcfb013df6c934ba3406b57..f69da67e2102f3855b1340518f6f5f39cbd2f053 100644
--- a/sdpm/src/sdpmSolverTransonic.h
+++ b/sdpm/src/sdpmSolverTransonic.h
@@ -35,7 +35,7 @@ class SDPM_API SolverTransonic : public Solver
sdpmVectorXd _dCorr; ///< transonic correction terms
public:
- SolverTransonic(Problem &pbl);
+ SolverTransonic(Problem &pbl, int vrb = 1);
void run() override;
diff --git a/sdpm/tests/lann_xsonic.py b/sdpm/tests/lann_xsonic.py
index 9f24e3438d9889343a166f1fb3c37f62e5fe0991..3480595cde561997a8fa75cf452d370d797d4b58 100644
--- a/sdpm/tests/lann_xsonic.py
+++ b/sdpm/tests/lann_xsonic.py
@@ -53,7 +53,7 @@ def main():
wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, 1, 1)
wing.addMotion()
wing.setRigidMotion(0, amp, 0., xf/math.sqrt(1-minf*minf), 0.)
- dcp = interpolate_pressure(build_fpath('models/lann_dcp_rans', __file__, 1), wing.getNodes())
+ dcp = interpolate_pressure(build_fpath('models/lann_dcp_rans.csv', __file__, 1), wing.getNodes())
wing.setTransonicPressureGrad(dcp[:,0])
pbl.addBody(wing)
tms['pbl'].stop()
diff --git a/sdpm/utils.py b/sdpm/utils.py
index 169b9eafa45b24201a90232c1df6e3e6162cd237..7e432f6c61759300204543d8cea6694eec55e24a 100644
--- a/sdpm/utils.py
+++ b/sdpm/utils.py
@@ -158,7 +158,7 @@ def interpolate_pressure(fname, surf_nodes):
Parameters
----------
fname : string
- name of the .csv file (without extension)
+ name of the CSV file containing the pressure derivative
surf_nodes : array
nodes of surface grid
@@ -170,7 +170,7 @@ def interpolate_pressure(fname, surf_nodes):
import numpy as np
import scipy.interpolate as spip
# Read pressure and remove duplicates
- cp = np.loadtxt(fname+'.csv', delimiter=',', skiprows=1)
+ cp = np.loadtxt(fname, delimiter=',', skiprows=1)
_, idx = np.unique(cp[:,:3], axis=0, return_index=True)
cp = cp[idx,:]
# Interpolate on surface mesh