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

Add transonic solver in API.

parent e95ef0a3
No related branches found
No related tags found
1 merge request!3Version 1.2
Pipeline #23410 passed
...@@ -50,7 +50,7 @@ def init_sdpm(cfg, use_ad=False): ...@@ -50,7 +50,7 @@ def init_sdpm(cfg, use_ad=False):
# Imports # Imports
import sdpm import sdpm
from sdpm.gmsh import MeshLoader from sdpm.gmsh import MeshLoader
from sdpm.utils import interpolate_modes from sdpm.utils import interpolate_modes, interpolate_pressure
import math import math
# Sanity checks # Sanity checks
...@@ -79,11 +79,14 @@ def init_sdpm(cfg, use_ad=False): ...@@ -79,11 +79,14 @@ def init_sdpm(cfg, use_ad=False):
x_modes, amp_modes = interpolate_modes(cfg['Modes'], _bdy.getNodes()) x_modes, amp_modes = interpolate_modes(cfg['Modes'], _bdy.getNodes())
for i in range(cfg['Num_modes']): for i in range(cfg['Num_modes']):
_bdy.setFlexibleMotion(i, 1 / amp_modes[i], x_modes[:, 3*i:3*i+3]) _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) _pbl.addBody(_bdy)
# Solver # Solver
vrb = cfg['Verb_lvl'] if 'Verb_lvl' in cfg else 1 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 _adj = sdpm.Adjoint(_sol) if use_ad else None
# Return # Return
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
using namespace sdpm; 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 // Build node row id
int i = 0; int i = 0;
...@@ -40,9 +40,10 @@ void SolverTransonic::run() ...@@ -40,9 +40,10 @@ void SolverTransonic::run()
computeCorrection(); computeCorrection();
std::cout << "done." << std::endl; std::cout << "done." << std::endl;
std::cout << "Running transonic SDPM solver...\n" std::cout << "Running transonic SDPM solver..." << std::flush;
<< "-- Mean flow" << std::endl; if (_vrb > 0)
std::cout << "computing sources... " << std::flush; std::cout << "\n-- Mean flow\n"
<< "computing sources... " << std::flush;
sdpmDouble beta = _pbl.getCompressibility(); sdpmDouble beta = _pbl.getCompressibility();
sdpmVector3d ui = _pbl.getDragDir(); sdpmVector3d ui = _pbl.getDragDir();
sdpmVectorXd etauxy(_nP); sdpmVectorXd etauxy(_nP);
...@@ -55,23 +56,27 @@ void SolverTransonic::run() ...@@ -55,23 +56,27 @@ void SolverTransonic::run()
etauz(_rows[e]) = ui[2] * e->getCompressibleNormal()[2]; etauz(_rows[e]) = ui[2] * e->getCompressibleNormal()[2];
} }
} }
std::cout << "done." << std::endl; if (_vrb > 0)
std::cout << "solving for doublets... " << std::flush; std::cout << "done.\n"
<< "solving for doublets... " << std::flush;
sdpmVectorXd rhsxy = _B0 * etauxy; sdpmVectorXd rhsxy = _B0 * etauxy;
sdpmVectorXd rhsz = _B0 * etauz; sdpmVectorXd rhsz = _B0 * etauz;
sdpmVectorXd emuxy = solve(_A0, rhsxy); sdpmVectorXd emuxy = solve(_A0, rhsxy);
sdpmVectorXd emuz = solve(_A0, rhsz); sdpmVectorXd emuz = solve(_A0, rhsz);
std::cout << "done." << std::endl; if (_vrb > 0)
std::cout << "computing flow and loads on bodies... " << std::flush; std::cout << "done.\n"
<< "computing flow and loads on bodies... " << std::flush;
post(etauxy + etauz, emuxy + emuz.cwiseProduct(_dCorr)); 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 imd = 0; imd < _pbl.getNumModes(); ++imd)
{ {
for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq) for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
{ {
std::cout << "-- Mode #" << imd << ", frequency #" << ifq << std::endl; if (_vrb > 0)
std::cout << "computing sources... " << std::flush; std::cout << "-- Mode #" << imd << ", frequency #" << ifq << "\n"
<< "computing sources... " << std::flush;
sdpmVectorXcd etau1(_nP); sdpmVectorXcd etau1(_nP);
sdpmVectorXcd etau1xy(_nP); sdpmVectorXcd etau1xy(_nP);
sdpmVectorXcd etau1z(_nP); sdpmVectorXcd etau1z(_nP);
...@@ -88,9 +93,10 @@ void SolverTransonic::run() ...@@ -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); 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) // Split complex set of equations into real set (twice the size)
std::cout << "solving for doublets... " << std::flush;
sdpmMatrixXd A(2 * _nP, 2 * _nP); sdpmMatrixXd A(2 * _nP, 2 * _nP);
A.block(0, 0, _nP, _nP) = _A[ifq].real(); A.block(0, 0, _nP, _nP) = _A[ifq].real();
A.block(_nP, _nP, _nP, _nP) = _A[ifq].real(); A.block(_nP, _nP, _nP, _nP) = _A[ifq].real();
...@@ -106,22 +112,29 @@ void SolverTransonic::run() ...@@ -106,22 +112,29 @@ void SolverTransonic::run()
sdpmVectorXd muz = solve(A, rhs1z); sdpmVectorXd muz = solve(A, rhs1z);
sdpmVectorXcd emu1xy = muxy.block(0, 0, _nP, 1) + sdpmComplex(0, 1) * muxy.block(_nP, 0, _nP, 1); 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); sdpmVectorXcd emu1z = muz.block(0, 0, _nP, 1) + sdpmComplex(0, 1) * muz.block(_nP, 0, _nP, 1);
std::cout << "done." << std::endl; if (_vrb > 0)
std::cout << "computing flow and loads on bodies... " << std::flush; std::cout << "done.\n"
<< "computing flow and loads on bodies... " << std::flush;
sdpmVectorXcd mu1c = emu1xy + _dCorr.cwiseProduct(emu1z); sdpmVectorXcd mu1c = emu1xy + _dCorr.cwiseProduct(emu1z);
post(imd, ifq, etau1, mu1c); post(imd, ifq, etau1, mu1c);
std::cout << "done." << std::endl; if (_vrb > 0)
std::cout << "done." << std::endl;
} }
} }
// GAF // 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) 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); 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; std::cout << std::endl;
} }
......
...@@ -35,7 +35,7 @@ class SDPM_API SolverTransonic : public Solver ...@@ -35,7 +35,7 @@ class SDPM_API SolverTransonic : public Solver
sdpmVectorXd _dCorr; ///< transonic correction terms sdpmVectorXd _dCorr; ///< transonic correction terms
public: public:
SolverTransonic(Problem &pbl); SolverTransonic(Problem &pbl, int vrb = 1);
void run() override; void run() override;
......
...@@ -53,7 +53,7 @@ def main(): ...@@ -53,7 +53,7 @@ def main():
wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, 1, 1) wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, 1, 1)
wing.addMotion() wing.addMotion()
wing.setRigidMotion(0, amp, 0., xf/math.sqrt(1-minf*minf), 0.) 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]) wing.setTransonicPressureGrad(dcp[:,0])
pbl.addBody(wing) pbl.addBody(wing)
tms['pbl'].stop() tms['pbl'].stop()
......
...@@ -158,7 +158,7 @@ def interpolate_pressure(fname, surf_nodes): ...@@ -158,7 +158,7 @@ def interpolate_pressure(fname, surf_nodes):
Parameters Parameters
---------- ----------
fname : string fname : string
name of the .csv file (without extension) name of the CSV file containing the pressure derivative
surf_nodes : array surf_nodes : array
nodes of surface grid nodes of surface grid
...@@ -170,7 +170,7 @@ def interpolate_pressure(fname, surf_nodes): ...@@ -170,7 +170,7 @@ def interpolate_pressure(fname, surf_nodes):
import numpy as np import numpy as np
import scipy.interpolate as spip import scipy.interpolate as spip
# Read pressure and remove duplicates # 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) _, idx = np.unique(cp[:,:3], axis=0, return_index=True)
cp = cp[idx,:] cp = cp[idx,:]
# Interpolate on surface mesh # Interpolate on surface mesh
......
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