diff --git a/flow/_src/floww.h b/flow/_src/floww.h index b593315c070fe2fec17d6ef57ea9a6de9313d0c4..297a990acc42d4fe0a490a953198d8c79d9765f0 100644 --- a/flow/_src/floww.h +++ b/flow/_src/floww.h @@ -18,11 +18,7 @@ #include "wAssign.h" #include "wBlowing.h" #include "wBoundary.h" -#include "wFace.h" #include "wFreestream.h" -#include "wF0El.h" -#include "wF0Ps.h" -#include "wF1El.h" #include "wKutta.h" #include "wKuttaElement.h" #include "wMedium.h" diff --git a/flow/_src/floww.i b/flow/_src/floww.i index 4f7387e49a2610f9f26a7c37672b7d1903fe5d6d..3e969f4799c26c62549499d49f29ebdcf5fe7012 100644 --- a/flow/_src/floww.i +++ b/flow/_src/floww.i @@ -51,7 +51,6 @@ threads="1" %shared_ptr(flow::Initial); %shared_ptr(flow::Dirichlet); %shared_ptr(flow::Freestream); -%shared_ptr(flow::Slip); %shared_ptr(flow::Wake); %shared_ptr(flow::Kutta); %shared_ptr(flow::Blowing); @@ -74,31 +73,6 @@ threads="1" %include "wWakeElement.h" %include "wKuttaElement.h" -%feature("director") flow::F0PsC; -%pythonappend flow::F0PsC::F0PsC "self.__disown__()" // only for directors -%feature("director") flow::F0PsPhiInf; -%pythonappend flow::F0PsPhiInf::F0PsPhiInf "self.__disown__()" // only for directors -%nodefault flow::F0Ps; // is pure virtual but swig doesn't know it. -%include "wF0Ps.h" - -%feature("director") flow::F0ElRho; -%pythonappend flow::F0ElRho::F0ElRho "self.__disown__()" // only for directors -%feature("director") flow::F0ElRhoL; -%pythonappend flow::F0ElRhoL::F0ElRhoL "self.__disown__()" // only for directors -%feature("director") flow::F0ElMach; -%pythonappend flow::F0ElMach::F0ElMach "self.__disown__()" // only for directors -%feature("director") flow::F0ElMachL; -%pythonappend flow::F0ElMachL::F0ElMachL "self.__disown__()" // only for directors -%feature("director") flow::F0ElCp; -%pythonappend flow::F0ElCp::F0ElCp "self.__disown__()" // only for directors -%feature("director") flow::F0ElCpL; -%pythonappend flow::F0ElCpL::F0ElCpL "self.__disown__()" // only for directors -%include "wF0El.h" - -%feature("director") flow::F1ElVi; -%pythonappend flow::F1ElVi::F1ElVi "self.__disown__()" // only for directors -%include "wF1El.h" - %immutable flow::Problem::msh; // read-only variable (no setter) %include "wProblem.h" diff --git a/flow/benchmark/onera.py b/flow/benchmark/onera.py index 85aa46ab31925a95f54c0a8a0a7cdacaa6859c0d..9edc59a548c86d4e203c4faaf1a3cfdf2b6212c0 100644 --- a/flow/benchmark/onera.py +++ b/flow/benchmark/onera.py @@ -65,7 +65,6 @@ def main(): # define flow variables alpha = 3.06*np.pi/180 M_inf = 0.839 - M_crit = 5 # Squared critical Mach number (above which density is modified) dim = 3 # define dimension @@ -80,7 +79,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') tms['pre'].stop() # solve problem diff --git a/flow/cases/coyote.py b/flow/cases/coyote.py index cb31118464ffc1a023e61e1f0a54e06e1d103fd7..6994b90a98679691bff9bf300ece7730ea60b44c 100644 --- a/flow/cases/coyote.py +++ b/flow/cases/coyote.py @@ -71,7 +71,6 @@ def getParam(): p['LS_tol'] = 1e-6 # Tolerance on linesearch residual p['Max_it_LS'] = 10 # Linesearch maximum number of iterations p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased - p['M_crit'] = 5. # Critical Mach number above which density is damped return p def main(): diff --git a/flow/cases/n0012.py b/flow/cases/n0012.py index a106ca95ddf9792bc66d2840437f6836d385f330..52234831583ada6a7bb0b899988ec49142faee5e 100644 --- a/flow/cases/n0012.py +++ b/flow/cases/n0012.py @@ -50,7 +50,6 @@ def getParam(): p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual p['Max_it'] = 20 # Solver maximum number of iterations p['Relaxation'] = 0.7 # Relaxation parameter for Picard solver - p['M_crit'] = 5. # Critical Mach number above which density is damped return p def main(): diff --git a/flow/cases/n0012_3.py b/flow/cases/n0012_3.py index cc2c03442293f5d4cf6b31f46d9f3f24df2bbe36..c0d087e4dbd76bf3a1fb8ce65fca2895067a34a8 100644 --- a/flow/cases/n0012_3.py +++ b/flow/cases/n0012_3.py @@ -58,7 +58,6 @@ def getParam(): p['LS_tol'] = 1e-6 # Tolerance on linesearch residual p['Max_it_LS'] = 10 # Linesearch maximum number of iterations p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased - p['M_crit'] = 5. # Critical Mach number above which density is damped return p def main(): diff --git a/flow/cases/n64A410.py b/flow/cases/n64A410.py index e406b3340f4f6f63be9e846fdd2b96fe1a0ad697..3b8fcaa7b792e9ac0e0818ea3c08418e44858818 100644 --- a/flow/cases/n64A410.py +++ b/flow/cases/n64A410.py @@ -52,7 +52,6 @@ def getParam(): p['LS_tol'] = 1e-6 # Tolerance on linesearch residual p['Max_it_LS'] = 10 # Linesearch maximum number of iterations p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased - p['M_crit'] = 5. # Critical Mach number above which density is damped return p def main(): diff --git a/flow/cases/rae2822.py b/flow/cases/rae2822.py index c5b9f0f206616ada6416870e9776dbcec59afb0f..a447e18d0828d85fbdf6ce78cdc12e2b36c0c820 100644 --- a/flow/cases/rae2822.py +++ b/flow/cases/rae2822.py @@ -55,7 +55,6 @@ def getParam(): p['LS_tol'] = 1e-6 # Tolerance on linesearch residual p['Max_it_LS'] = 10 # Linesearch maximum number of iterations p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased - p['M_crit'] = 5. # Critical Mach number above which density is damped return p def main(): diff --git a/flow/cases/wbht.py b/flow/cases/wbht.py index 69b2c484f6213ec02b8294814feaaee26f9eb4ca..2ab10ec0b533675dec8169c0c734c119a2b0dbed 100644 --- a/flow/cases/wbht.py +++ b/flow/cases/wbht.py @@ -83,7 +83,6 @@ def getParam(): p['LS_tol'] = 1e-6 # Tolerance on linesearch residual p['Max_it_LS'] = 10 # Linesearch maximum number of iterations p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased - p['M_crit'] = 5. # Critical Mach number above which density is damped return p def main(): diff --git a/flow/cases/wht.py b/flow/cases/wht.py index 1efcd6505ad112f458b39806a48eb602cdf99c51..72deaf3e25f4080fc4575bca4eaba94422ee07b0 100644 --- a/flow/cases/wht.py +++ b/flow/cases/wht.py @@ -74,7 +74,6 @@ def getParam(): p['LS_tol'] = 1e-6 # Tolerance on linesearch residual p['Max_it_LS'] = 10 # Linesearch maximum number of iterations p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased - p['M_crit'] = 5. # Critical Mach number above which density is damped return p def main(): diff --git a/flow/default.py b/flow/default.py index a6b3e9b4e6591ca3c9b365a761dc17150530c1cb..50b4bb56069f773d187c19d7f2126cad2b001338 100644 --- a/flow/default.py +++ b/flow/default.py @@ -42,29 +42,26 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None): gmshWriter.save(msh.name) return msh, gmshWriter -def problem(msh, dim, alpha, beta, minf, mcrit, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = None, tp = None, vsc = False, dbc = False): +def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = None, tp = None, vsc = False, dbc = False): '''Initialize problem ''' # create base formulation pbl = flo.Problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref) # add incompressible or compressible air medium - if minf == 0: - pbl.set(flo.Medium(msh, fld, flo.F0ElRhoL(), flo.F0ElMachL(), flo.F0ElCpL(), flo.F0PsPhiInf(dim, alpha, beta))) - else: - pbl.set(flo.Medium(msh, fld, flo.F0ElRho(minf, mcrit), flo.F0ElMach(minf), flo.F0ElCp(minf), flo.F0PsPhiInf(dim, alpha, beta))) + pbl.set(flo.Medium(msh, fld, minf, dim, alpha, beta)) # add initial condition - pbl.set(flo.Initial(msh, fld, flo.F0PsPhiInf(dim, alpha, beta))) + pbl.set(flo.Initial(msh, fld, dim, alpha, beta)) if dbc: # Dirichlet BCs provided and explicitely imposed - dirichlet = flo.Dirichlet(msh, upstr, flo.F0PsPhiInf(dim, alpha, beta)) + dirichlet = flo.Dirichlet(msh, upstr, dim, alpha, beta) pbl.add(dirichlet) else: # Dirichlet BCs not provided, Neumann imposed instead, and a DOF will be pinned automatically dirichlet = None - pbl.add(flo.Freestream(msh, upstr, flo.F1ElVi(dim, alpha, beta))) + pbl.add(flo.Freestream(msh, upstr, dim, alpha, beta)) # add freestream (Neumann) conditions - pbl.add(flo.Freestream(msh, ff, flo.F1ElVi(dim, alpha, beta))) - pbl.add(flo.Freestream(msh, dnstr, flo.F1ElVi(dim, alpha, beta))) + pbl.add(flo.Freestream(msh, ff, dim, alpha, beta)) + pbl.add(flo.Freestream(msh, dnstr, dim, alpha, beta)) # add wake and Kutta conditions if dim == 2 and te is not None: wake = flo.Wake(msh, [wk, wk+'_', fld]) diff --git a/flow/scripts/config.py b/flow/scripts/config.py index 8625c7a4591c1d91fb0c63dbb484b9037a782a16..93c79c755c5c3f468a910f40ae1ffeaf0f38e9fa 100644 --- a/flow/scripts/config.py +++ b/flow/scripts/config.py @@ -19,6 +19,7 @@ ## Base class for configuration # Adrien Crovato +import math import flow as f import tbox import tbox.gmsh as gmsh @@ -39,11 +40,16 @@ class Config: if p['Dim'] == 2: if 'AoS' in p and p['AoS'] != 0: raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n') + beta = 0 if 'Symmetry' in p: raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n') else: - if 'AoS' in p and p['AoS'] != 0 and 'Symmetry' in p: - raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n') + if 'AoS' in p: + if p['AoS'] != 0 and 'Symmetry' in p: + raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n') + beta = p['AoS']*math.pi/180 + else: + beta = 0 # basic config if p['Format'] == 'vtk': try: @@ -82,19 +88,14 @@ class Config: print('\n') # initialize the problem - self.pbl = f.Problem(self.msh, p['Dim'], 0., 0., p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref']) - self.phiInfFun = f.F0PsPhiInf(p['Dim'], 0.) - self.velInfFun = f.F1ElVi(p['Dim'], 0.) + self.pbl = f.Problem(self.msh, p['Dim'], 0., beta, p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref']) # add a medium "air" - if p['M_inf'] == 0: - self.pbl.set(f.Medium(self.msh, p['Fluid'], f.F0ElRhoL(), f.F0ElMachL(), f.F0ElCpL(), self.phiInfFun)) - else: - self.pbl.set(f.Medium(self.msh, p['Fluid'], f.F0ElRho(p['M_inf'], p['M_crit']), f.F0ElMach(p['M_inf']), f.F0ElCp(p['M_inf']), self.phiInfFun)) + self.pbl.set(f.Medium(self.msh, p['Fluid'], p['M_inf'], p['Dim'], 0., beta)) # add initial condition - self.pbl.set(f.Initial(self.msh, p['Fluid'], self.phiInfFun)) + self.pbl.set(f.Initial(self.msh, p['Fluid'], p['Dim'], 0., beta)) # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically) for i in range (0, len(p['Farfield'])): - self.pbl.add(f.Freestream(self.msh, p['Farfield'][i], self.velInfFun)) + self.pbl.add(f.Freestream(self.msh, p['Farfield'][i], p['Dim'], 0., beta)) # add solid boundaries if p['Dim'] == 2: self.bnd = f.Boundary(self.msh, [p['Wing'], p['Fluid']]) diff --git a/flow/scripts/polar.py b/flow/scripts/polar.py index 7cb623ee92de2093affcac9941ff4c4bc392e9c2..8d3e694275988ebf54676c3dc5a2543e9a5a241f 100644 --- a/flow/scripts/polar.py +++ b/flow/scripts/polar.py @@ -40,10 +40,6 @@ class Polar(Config): self.tag = None self.mach = p['M_inf'] self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step']) - if 'AoS' in p: - self.beta = p['AoS']*math.pi/180 - else: - self.beta = 0 # basic init Config.__init__(self, p) @@ -58,10 +54,7 @@ class Polar(Config): # define current angle of attack alpha = alpha*math.pi/180 # update problem - self.pbl.alpha = alpha - self.pbl.beta = self.beta - self.phiInfFun.update(alpha, self.beta) - self.velInfFun.update(alpha, self.beta) + self.pbl.update(alpha) # run the solver and save the results print(ccolors.ANSI_BLUE + 'Running the solver for', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET) self.tms["solver"].start() diff --git a/flow/scripts/trim.py b/flow/scripts/trim.py index a14bea0cf58031af1aeb04727bee201c6d2a5915..1c3b0f5f55f7a3b4122b07507b696d4b35b89ca2 100644 --- a/flow/scripts/trim.py +++ b/flow/scripts/trim.py @@ -41,10 +41,6 @@ class Trim(Config): self.mach = p['M_inf'] self.ClTgt = p['CL'] self.alpha = p['AoA']*math.pi/180 - if 'AoS' in p: - self.beta = p['AoS']*math.pi/180 - else: - self.beta = 0 self.dCl = p['dCL'] # basic init Config.__init__(self, p) @@ -56,10 +52,7 @@ class Trim(Config): # define current angle of attack print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET) # update problem - self.pbl.alpha = self.alpha - self.pbl.beta = self.beta - self.phiInfFun.update(self.alpha, self.beta) - self.velInfFun.update(self.alpha, self.beta) + self.pbl.update(self.alpha) # run the solver self.tms["solver"].start() success = self.solver.run() diff --git a/flow/src/flow.h b/flow/src/flow.h index 351ca549ddaf0367b01d74162f0a773b07b5060f..27277a2b0274dd6320928bb558e633bb924a640a 100644 --- a/flow/src/flow.h +++ b/flow/src/flow.h @@ -36,7 +36,6 @@ */ namespace flow { - // Problem and B.C. handling class Problem; class Medium; @@ -69,7 +68,19 @@ class Adjoint; // General class F0Ps; class F0El; +class F0ElC; +class F0ElRho; +class F0ElRhoL; +class F0ElMach; +class F0ElMachL; +class F0ElCp; +class F0ElCpL; class F1El; +class F1ElVi; +class F1Ct; +class F1CtDrag; +class F1CtSide; +class F1CtLift; class Face; class FaceEq; diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp index 11eebf9865a1b6cad58ff8ec79553f6f884c2c13..3eac871cce555e9f15a71f6f84af641dbe964367 100644 --- a/flow/src/wAdjoint.cpp +++ b/flow/src/wAdjoint.cpp @@ -247,28 +247,22 @@ void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD) // Multithread tbb::spin_mutex mutex; - // Reset gradients + // Reset and compute gradient of lift and drag dL.setZero(); dD.setZero(); - - // Lift and drag normalized directions - Eigen::RowVector3d dirL, dirD; - std::tie(dirD, std::ignore, dirL) = sol->pbl->computeDirections(); - - // Gradient of lift and drag for (auto sur : sol->pbl->bnds) { tbb::parallel_do(sur->groups[0]->tag->elems.begin(), sur->groups[0]->tag->elems.end(), [&](Element *e) { // Associated volume element Element *eV = sur->svMap.at(e); // Build - Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, sol->pbl->medium->cP); + Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, *sol->pbl->medium->cP); // Assembly tbb::spin_mutex::scoped_lock lock(mutex); for (size_t i = 0; i < e->nodes.size(); ++i) { - Eigen::RowVectorXd AeL = dirL * Ae.block(i * 3, 0, 3, eV->nodes.size()); - Eigen::RowVectorXd AeD = dirD * Ae.block(i * 3, 0, 3, eV->nodes.size()); + Eigen::RowVectorXd AeL = sol->pbl->dirL.eval().transpose() * Ae.block(i * 3, 0, 3, eV->nodes.size()); + Eigen::RowVectorXd AeD = sol->pbl->dirD.eval().transpose() * Ae.block(i * 3, 0, 3, eV->nodes.size()); for (size_t j = 0; j < eV->nodes.size(); ++j) { Node *nodj = eV->nodes[j]; @@ -307,26 +301,6 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR) } }); } - - // Wake contribution (through stabilisation term) - // @todo very slight influence, can be removed by using Vi = [1,0,0]. Same results, same cost, less general but less code - for (auto wake : sol->pbl->wBCs) - { - Eigen::VectorXd dvi = sol->pbl->computeGradientVi(); // gradient of freestream velocity - tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) { - // Build bu, bl - Eigen::VectorXd bue, ble; - std::tie(bue, ble) = WakeResidual::build(*we, sol->phi, dvi); - // Assembly - tbb::spin_mutex::scoped_lock lock(mutex); - for (size_t ii = 0; ii < we->nRow; ++ii) - { - Node *nodi = we->wkN[ii].second; - dR(nodi->row) += bue(ii); - dR(nodi->row) -= ble(ii); - } - }); - } } /** @@ -335,17 +309,14 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR) */ void Adjoint::buildGradientLoadsAoa(double &dL, double &dD) { - // Gradient of lift and drag normalized directions - Eigen::RowVector3d dirL, dirD; - std::tie(dirD, std::ignore, dirL) = sol->pbl->computeGradientDirections(); // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area) Eigen::Vector3d Cf(0, 0, 0); for (auto bnd : sol->pbl->bnds) for (size_t i = 0; i < bnd->nodes.size(); ++i) Cf += bnd->nLoads[i]; // Rotate to gradient of flow direction - dD = Cf.dot(dirD); - dL = Cf.dot(dirL); + dD = Cf.dot(sol->pbl->dirD.evalGrad()); + dL = Cf.dot(sol->pbl->dirL.evalGrad()); } /** @@ -388,7 +359,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa } } // Assembly (supersonic) - double mach = fluid->mach.eval(*e, sol->phi, 0); + double mach = fluid->mach->eval(*e, sol->phi, 0); if (mach > sol->mCOv) { for (size_t i = 0; i < e->nodes.size(); ++i) @@ -413,11 +384,10 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa // Wake for (auto wake : sol->pbl->wBCs) { - Eigen::VectorXd vi = sol->pbl->computeVi(); tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) { // Build elementary matrices Eigen::MatrixXd Aeu, Ael, Aes; - std::tie(Aeu, Ael, Aes) = WakeResidual::buildGradientMesh(*we, sol->phi, vi); + std::tie(Aeu, Ael, Aes) = WakeResidual::buildGradientMesh(*we, sol->phi); // Assembly tbb::spin_mutex::scoped_lock lock(mutex); for (size_t i = 0; i < we->nRow; ++i) @@ -511,15 +481,9 @@ void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd // Multithread tbb::spin_mutex mutex; - // Reset gradients + // Reset and compute gradient of lift and drag dL.setZero(); dD.setZero(); - - // Lift and drag normalized directions - Eigen::RowVector3d dirL, dirD; - std::tie(dirD, std::ignore, dirL) = sol->pbl->computeDirections(); - - // Gradient of lift and drag for (auto bnd : sol->pbl->bnds) { tbb::parallel_do(bnd->groups[0]->tag->elems.begin(), bnd->groups[0]->tag->elems.end(), [&](Element *e) { @@ -527,16 +491,16 @@ void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd Element *eV = bnd->svMap.at(e); // Build Eigen::MatrixXd Aes, Aev; - std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, sol->pbl->medium->cP, sol->pbl->nDim); + std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, *sol->pbl->medium->cP, sol->pbl->nDim); // Assembly tbb::spin_mutex::scoped_lock lock(mutex); for (size_t i = 0; i < e->nodes.size(); ++i) { // rotate to flow direction - Eigen::RowVectorXd AesL = dirL * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size()); - Eigen::RowVectorXd AesD = dirD * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size()); - Eigen::RowVectorXd AevL = dirL * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size()); - Eigen::RowVectorXd AevD = dirD * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size()); + Eigen::RowVectorXd AesL = sol->pbl->dirL.eval().transpose() * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size()); + Eigen::RowVectorXd AesD = sol->pbl->dirD.eval().transpose() * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size()); + Eigen::RowVectorXd AevL = sol->pbl->dirL.eval().transpose() * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size()); + Eigen::RowVectorXd AevD = sol->pbl->dirD.eval().transpose() * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size()); for (size_t m = 0; m < sol->pbl->nDim; ++m) { // surface diff --git a/flow/src/wAssign.cpp b/flow/src/wAssign.cpp index 7056fcb68be743a1ed80ba117571d2a8fd640515..82ecc405ec31ee9d38e93a5ccd24833cb132e5b9 100644 --- a/flow/src/wAssign.cpp +++ b/flow/src/wAssign.cpp @@ -23,30 +23,31 @@ using namespace tbox; using namespace flow; -Assign::Assign(std::shared_ptr<MshData> _msh, int no, F0Ps &_f) : Group(_msh, no), f(_f) +Assign::Assign(std::shared_ptr<MshData> _msh, int no) : Group(_msh, no) { getNodes(); } -Assign::Assign(std::shared_ptr<MshData> _msh, std::string const &name, F0Ps &_f) : Group(_msh, name), f(_f) +Assign::Assign(std::shared_ptr<MshData> _msh, std::string const &name) : Group(_msh, name) { getNodes(); } +Assign::~Assign() +{ + delete f; + std::cout << "~Assign()\n"; +} /** * @brief Prescribe value to all nodes (default implementation) - * - * Value can be constant or position dependent - * @authors Adrien Crovato, Romain Boman */ void Assign::apply(std::vector<double> &vec) { for (auto n : nodes) - vec[n->row] = f.eval(n->pos); + vec[n->row] = f->eval(n->pos); } /** * @brief Get all nodes belonging to boundary tag - * @authors Adrien Crovato, Romain Boman */ void Assign::getNodes() { @@ -58,31 +59,35 @@ void Assign::getNodes() nodes.resize(std::distance(nodes.begin(), it)); } -Initial::Initial(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f) : Assign(_msh, no, _f) +Initial::Initial(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta) : Assign(_msh, no) { + f = new F0PsPhiInf(dim, alpha, beta); } -Initial::Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f) : Assign(_msh, name, _f) +Initial::Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Assign(_msh, name) { + f = new F0PsPhiInf(dim, alpha, beta); } void Initial::write(std::ostream &out) const { out << "flow::Initial condition on " << *tag << ")\n"; for (auto n : nodes) - out << '\t' << *n << ": (val=" << f.eval(n->pos) << ")\n"; + out << '\t' << *n << ": (val=" << f->eval(n->pos) << ")\n"; } -Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f, bool pin) : Assign(_msh, no, _f) +Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta, bool pin) : Assign(_msh, no) { + f = new F0PsPhiInf(dim, alpha, beta); // Only retain the first node, so that the DOF associated to this node will be pinned if (pin) this->nodes.resize(1); } -Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f) : Assign(_msh, name, _f) +Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Assign(_msh, name) { + f = new F0PsPhiInf(dim, alpha, beta); } void Dirichlet::write(std::ostream &out) const { out << "flow::Dirichlet boundary condition on " << *tag << ")\n"; for (auto n : nodes) - out << '\t' << *n << ": (val=" << f.eval(n->pos) << ")\n"; + out << '\t' << *n << ": (val=" << f->eval(n->pos) << ")\n"; } \ No newline at end of file diff --git a/flow/src/wAssign.h b/flow/src/wAssign.h index e86d793b40d1da0b23e2e50b1f650cc023c544e5..784c51e9d37e1adbdb83b3b36faba62bf632eb58 100644 --- a/flow/src/wAssign.h +++ b/flow/src/wAssign.h @@ -36,11 +36,11 @@ class FLOW_API Assign : public tbox::Group public: std::vector<tbox::Node *> nodes; ///< nodes of the boundary #ifndef SWIG - F0Ps &f; ///< position-based function to compute the boundary value + F0Ps *f; ///< position-based function to compute the boundary value - Assign(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f); - Assign(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f); - virtual ~Assign() { std::cout << "~Assign()\n"; } + Assign(std::shared_ptr<tbox::MshData> _msh, int no); + Assign(std::shared_ptr<tbox::MshData> _msh, std::string const &name); + virtual ~Assign(); void apply(std::vector<double> &vec); #endif @@ -56,8 +56,8 @@ private: class FLOW_API Initial : public Assign { public: - Initial(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f); - Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f); + Initial(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta = 0.0); + Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta = 0.0); virtual ~Initial() { std::cout << "~Initial()\n"; } #ifndef SWIG @@ -72,8 +72,8 @@ public: class FLOW_API Dirichlet : public Assign { public: - Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f, bool pin = false); - Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f); + Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta = 0.0, bool pin = false); + Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta = 0.0); virtual ~Dirichlet() { std::cout << "~Dirichlet()\n"; } #ifndef SWIG diff --git a/flow/src/wBlowing.cpp b/flow/src/wBlowing.cpp index 29420b69af6ad3fce9274cb3d17f8a4eed4cdc3b..509924621855ff2a78abe01701b40099339f0d22 100644 --- a/flow/src/wBlowing.cpp +++ b/flow/src/wBlowing.cpp @@ -39,7 +39,6 @@ Blowing::~Blowing() /** * @brief Create data structure - * @authors Adrien Crovato */ void Blowing::create() { @@ -55,14 +54,13 @@ void Blowing::create() uE.resize(tag->elems.size()); for (size_t i = 0; i < tag->elems.size(); ++i) { - Fct0C *f = new Fct0C(0.); - uE[i] = std::pair<Element *, Fct0C *>(tag->elems[i], f); + F0ElC *f = new F0ElC(0.); + uE[i] = std::pair<Element *, F0ElC *>(tag->elems[i], f); } } /** * @brief Set blowing velocity - * @authors Adrien Crovato */ void Blowing::setU(size_t i, double ue) { diff --git a/flow/src/wBlowing.h b/flow/src/wBlowing.h index de0f64007ef14a45b1e12cbb1d6a7bacd6141dee..c3706c31888a80a124d2500a707a38d17f680b00 100644 --- a/flow/src/wBlowing.h +++ b/flow/src/wBlowing.h @@ -19,7 +19,7 @@ #include "flow.h" #include "wGroup.h" -#include "wFct0.h" +#include "wF0El.h" #include <vector> namespace flow @@ -34,7 +34,7 @@ class FLOW_API Blowing : public tbox::Group public: std::vector<tbox::Node *> nodes; ///< nodes of Blowing #ifndef SWIG - std::vector<std::pair<tbox::Element *, tbox::Fct0C *>> uE; ///< blowing velocity for each element + std::vector<std::pair<tbox::Element *, F0ElC *>> uE; ///< blowing velocity for each element #endif Blowing(std::shared_ptr<tbox::MshData> _msh, int nos); diff --git a/flow/src/wBlowingResidual.cpp b/flow/src/wBlowingResidual.cpp index 150924687d7b1185ae5fc47dfa36c180a70482cb..7753c35891c498f8c1368d26648fbd2a5a24f6c3 100644 --- a/flow/src/wBlowingResidual.cpp +++ b/flow/src/wBlowingResidual.cpp @@ -15,7 +15,7 @@ */ #include "wBlowingResidual.h" -#include "wFct0.h" +#include "wF0El.h" #include "wElement.h" #include "wCache.h" @@ -29,7 +29,7 @@ using namespace flow; * * b = \int psi * vn dV */ -Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, Fct0 const &f) +Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, F0El const &f) { // Get pre-computed values Cache &cache = e.getVCache(); diff --git a/flow/src/wBlowingResidual.h b/flow/src/wBlowingResidual.h index 95e09a3e1b12f7ade26c923324c344fed722a424..8c77ef557a8470d9227e1c5b86dae57ed8ca8a0b 100644 --- a/flow/src/wBlowingResidual.h +++ b/flow/src/wBlowingResidual.h @@ -32,7 +32,7 @@ class FLOW_API BlowingResidual { public: // Newton - static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, tbox::Fct0 const &f); + static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, F0El const &f); }; } // namespace flow diff --git a/flow/src/wF0El.cpp b/flow/src/wF0El.cpp index 5643219538e53a449da43f1c0d6ca8f020c465dc..0829796ea3c936f49a438b3d8b85b847e61a7652 100644 --- a/flow/src/wF0El.cpp +++ b/flow/src/wF0El.cpp @@ -18,9 +18,30 @@ using namespace tbox; using namespace flow; +/** + * @brief Update the value + */ +void F0ElC::update(double _v) +{ + v = _v; +} +/** + * @brief Evaluate the constant + */ +double F0ElC::eval(Element const &e, std::vector<double> const &u, size_t k) const +{ + return v; +} +/** + * @brief Evaluate the gradient of the constant + */ +double F0ElC::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const +{ + return 0; +} + /** * @brief Evaluate the nonlinear density (constant over element) - * @authors Adrien Crovato */ double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) const { @@ -41,9 +62,8 @@ double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) c } /** * @brief Evaluate the nonlinear density derivative (constant over element) - * @authors Adrien Crovato */ -double F0ElRho::evalD(Element const &e, std::vector<double> const &u, size_t k) const +double F0ElRho::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const { double gradU = e.computeGradient(u, k).norm(); // norm of velocity @@ -60,14 +80,9 @@ double F0ElRho::evalD(Element const &e, std::vector<double> const &u, size_t k) return -sqrt(a * a * a * a * a) / ((1 + beta * (gradU / gradUC - 1)) * (1 + beta * (gradU / gradUC - 1))) * beta / gradUC / gradU; // has to be multiplied by grad u to recover true derivative } } -void F0ElRho::write(std::ostream &out) const -{ - out << "F0ElRho"; -} /** * @brief Evaluate the linear density (constant over element) - * @authors Adrien Crovato */ double F0ElRhoL::eval(Element const &e, std::vector<double> const &u, size_t k) const { @@ -75,20 +90,14 @@ double F0ElRhoL::eval(Element const &e, std::vector<double> const &u, size_t k) } /** * @brief Evaluate the linear density derivative (constant over element) - * @authors Adrien Crovato */ -double F0ElRhoL::evalD(Element const &e, std::vector<double> const &u, size_t k) const +double F0ElRhoL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const { return 0.; } -void F0ElRhoL::write(std::ostream &out) const -{ - out << "F0ElRhoL"; -} /** * @brief Evaluate the nonlinear Mach number (constant over element) - * @authors Adrien Crovato */ double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k) const { @@ -100,9 +109,8 @@ double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k) } /** * @brief Evaluate the nonlinear Mach number derivative (constant over element) - * @authors Adrien Crovato */ -double F0ElMach::evalD(Element const &e, std::vector<double> const &u, size_t k) const +double F0ElMach::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const { Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity double gradU2 = gradU.squaredNorm(); // velocity squared @@ -110,14 +118,9 @@ double F0ElMach::evalD(Element const &e, std::vector<double> const &u, size_t k) return 1 / sqrt(gradU2 * a2) + 0.5 * (gamma - 1) * sqrt(gradU2) / sqrt(a2 * a2 * a2); // has to be multiplied by grad u to recover true derivative } -void F0ElMach::write(std::ostream &out) const -{ - out << "F0ElMach"; -} /** * @brief Evaluate the linear Mach number (constant over element) - * @authors Adrien Crovato */ double F0ElMachL::eval(Element const &e, std::vector<double> const &u, size_t k) const { @@ -125,20 +128,14 @@ double F0ElMachL::eval(Element const &e, std::vector<double> const &u, size_t k) } /** * @brief Evaluate the linear Mach number derivative (constant over element) - * @authors Adrien Crovato */ -double F0ElMachL::evalD(Element const &e, std::vector<double> const &u, size_t k) const +double F0ElMachL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const { return 0.; } -void F0ElMachL::write(std::ostream &out) const -{ - out << "F0ElMachL"; -} /** * @brief Evaluate the nonlinear pressure coefficient (constant over element) - * @authors Adrien Crovato */ double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) const { @@ -154,9 +151,8 @@ double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) co } /** * @brief Evaluate the nonlinear pressure coefficient derivative (constant over element) - * @authors Adrien Crovato */ -double F0ElCp::evalD(Element const &e, std::vector<double> const &u, size_t k) const +double F0ElCp::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const { Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity double gradU2 = gradU.squaredNorm(); // velocity squared @@ -168,14 +164,9 @@ double F0ElCp::evalD(Element const &e, std::vector<double> const &u, size_t k) c else return -2 * sqrt(a * a * a * a * a); // has to be multiplied by grad u to recover true derivative } -void F0ElCp::write(std::ostream &out) const -{ - out << "F0ElCp"; -} /** * @brief Evaluate the linear pressure coefficient (constant over element) - * @authors Adrien Crovato */ double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) const { @@ -183,11 +174,7 @@ double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) c double gradU2 = gradU.squaredNorm(); // velocity squared return 1 - gradU2; } -double F0ElCpL::evalD(Element const &e, std::vector<double> const &u, size_t k) const +double F0ElCpL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const { return -2.; // has to be multiplied by grad u to recover true derivative -} -void F0ElCpL::write(std::ostream &out) const -{ - out << "F0ElCpL"; } \ No newline at end of file diff --git a/flow/src/wF0El.h b/flow/src/wF0El.h index a92e4d4b5410eb3eea924ac8bcfec9f2987ade06..d046d9919a26aca1055c69564eba8e2a28eb5bc1 100644 --- a/flow/src/wF0El.h +++ b/flow/src/wF0El.h @@ -18,153 +18,149 @@ #define WF0EL_H #include "flow.h" -#include "wFct0.h" +#include "wObserver.h" #include "wElement.h" #include <vector> namespace flow { +/** + * @brief Scalar function to be integrated over an element + */ +class FLOW_API F0El : public fwk::Observer +{ +public: + F0El() : Observer() {} + virtual ~F0El() {} + + virtual double eval(tbox::Element const &e, std::vector<double> const &u, + size_t k) const = 0; + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const = 0; +}; + +/** + * @brief Constant function + */ +class FLOW_API F0ElC : public F0El +{ + double v; + +public: + F0ElC(double _v) : F0El(), v(_v) {} + virtual void update(double _v) override; + + virtual double eval(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; +}; + /** * @brief Nonlinear density - * + * * The density is limited (Padé approximation) for large value of velocity and * particularized for diatomic gas. The input _mC2 is the square of the critical * Mach number used to compute the critical velocity. - * @authors Adrien Crovato */ -class FLOW_API F0ElRho : public tbox::Fct0 +class FLOW_API F0ElRho : public F0El { double gamma; ///< heat capacity ratio (diatomic gas only) double mInf; ///< freestream Mach number double gradUC; ///< critical velocity public: - F0ElRho(double _mInf, double _mC2) : Fct0(), gamma(1.4), mInf(_mInf) + F0ElRho(double _mInf, double _mC2 = 5) : F0El(), gamma(1.4), mInf(_mInf) { gradUC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2)); } -#ifndef SWIG - virtual double eval(tbox::Element const &e, - std::vector<double> const &u, + virtual double eval(tbox::Element const &e, std::vector<double> const &u, size_t k) const override; - virtual double evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; }; /** * @brief Linear density (constant) - * @authors Adrien Crovato */ -class FLOW_API F0ElRhoL : public tbox::Fct0 +class FLOW_API F0ElRhoL : public F0El { public: - F0ElRhoL() : Fct0() {} + F0ElRhoL() : F0El() {} -#ifndef SWIG - virtual double eval(tbox::Element const &e, - std::vector<double> const &u, + virtual double eval(tbox::Element const &e, std::vector<double> const &u, size_t k) const override; - virtual double evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; }; /** * @brief Nonlinear Mach number - * + * * Particularized for diatomic gas - * @authors Adrien Crovato */ -class FLOW_API F0ElMach : public tbox::Fct0 +class FLOW_API F0ElMach : public F0El { double gamma; ///< heat capacity ratio (diatomic gas only) double mInf; ///< freestream Mach number public: - F0ElMach(double _mInf) : Fct0(), gamma(1.4), mInf(_mInf) {} + F0ElMach(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {} -#ifndef SWIG - virtual double eval(tbox::Element const &e, - std::vector<double> const &u, + virtual double eval(tbox::Element const &e, std::vector<double> const &u, size_t k) const override; - virtual double evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; }; /** * @brief Linear Mach number (constant) - * @authors Adrien Crovato */ -class FLOW_API F0ElMachL : public tbox::Fct0 +class FLOW_API F0ElMachL : public F0El { public: - F0ElMachL() : Fct0() {} + F0ElMachL() : F0El() {} -#ifndef SWIG - virtual double eval(tbox::Element const &e, - std::vector<double> const &u, + virtual double eval(tbox::Element const &e, std::vector<double> const &u, size_t k) const override; - virtual double evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; }; /** * @brief Nonlinear pressure coefficient - * + * * Particularized for diatomic gas - * @authors Adrien Crovato */ -class FLOW_API F0ElCp : public tbox::Fct0 +class FLOW_API F0ElCp : public F0El { double gamma; ///< heat capacity ratio (diatomic gas only) double mInf; ///< freestream Mach number public: - F0ElCp(double _mInf) : Fct0(), gamma(1.4), mInf(_mInf) {} + F0ElCp(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {} -#ifndef SWIG - virtual double eval(tbox::Element const &e, - std::vector<double> const &u, + virtual double eval(tbox::Element const &e, std::vector<double> const &u, size_t k) const override; - virtual double evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; }; /** * @brief Linear pressure coefficient - * @authors Adrien Crovato */ -class FLOW_API F0ElCpL : public tbox::Fct0 +class FLOW_API F0ElCpL : public F0El { public: - F0ElCpL() : Fct0() {} + F0ElCpL() : F0El() {} -#ifndef SWIG - virtual double eval(tbox::Element const &e, - std::vector<double> const &u, + virtual double eval(tbox::Element const &e, std::vector<double> const &u, size_t k) const override; - virtual double evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; }; } // namespace flow diff --git a/flow/src/wF0Ps.cpp b/flow/src/wF0Ps.cpp index bc25ca8c4b4635a361772249e76d205e9f4b094f..7c8a0425d366ea94651cd3585a072bbbcfa11917 100644 --- a/flow/src/wF0Ps.cpp +++ b/flow/src/wF0Ps.cpp @@ -15,23 +15,9 @@ */ #include "wF0Ps.h" -#include <sstream> using namespace tbox; using namespace flow; -/** - * @brief Return constant function - * @authors Adrien Crovato - */ -double F0PsC::eval(Eigen::Vector3d const &pos) const -{ - return v; -} -void F0PsC::write(std::ostream &out) const -{ - out << "F0PsC"; -} - F0PsPhiInf::F0PsPhiInf(int _nDim, double _alpha, double _beta) : F0Ps(), nDim(_nDim), alpha(_alpha), @@ -48,17 +34,14 @@ F0PsPhiInf::F0PsPhiInf(int _nDim, double _alpha, double _beta) : F0Ps(), throw std::runtime_error("F0PsPhiInf::beta can be nonzero only for 3D cases!\n"); } /** - * @brief Update the angle of attack and the angle of sideslip - * @authors Adrien Crovato + * @brief Update the angle of attack */ -void F0PsPhiInf::update(double _alpha, double _beta) +void F0PsPhiInf::update(double _alpha) { alpha = _alpha; - beta = _beta; } /** * @brief Evaluate the freestream potential - * @authors Adrien Crovato */ double F0PsPhiInf::eval(Eigen::Vector3d const &pos) const { @@ -67,7 +50,3 @@ double F0PsPhiInf::eval(Eigen::Vector3d const &pos) const else return cos(alpha) * cos(beta) * pos(0) + sin(beta) * pos(1) + sin(alpha) * cos(beta) * pos(2); } -void F0PsPhiInf::write(std::ostream &out) const -{ - out << "F0PsPhiInf"; -} diff --git a/flow/src/wF0Ps.h b/flow/src/wF0Ps.h index b8b176075821f169bc7996dba989b17f010ddad9..b68f2dd8bbf8683fc6a1b469ab20610a4f56f04e 100644 --- a/flow/src/wF0Ps.h +++ b/flow/src/wF0Ps.h @@ -18,7 +18,7 @@ #define WF0PS_H #include "flow.h" -#include "wObject.h" +#include "wObserver.h" #include <Eigen/Dense> namespace flow @@ -26,44 +26,18 @@ namespace flow /** * @brief Scalar functions depending on spatial position - * @authors Adrien Crovato */ -class FLOW_API F0Ps : public fwk::wObject +class FLOW_API F0Ps : public fwk::Observer { public: -#ifndef SWIG - F0Ps() : wObject() - { - } -#endif - virtual ~F0Ps() - { - } -#ifndef SWIG - virtual double eval(Eigen::Vector3d const &pos) const = 0; -#endif -}; - -/** - * @brief Constant function - * @authors Adrien Crovato - */ -class FLOW_API F0PsC : public F0Ps -{ - double v; - -public: - F0PsC(double _v) : F0Ps(), v(_v) {} + F0Ps() : Observer() {} + virtual ~F0Ps() {} -#ifndef SWIG - virtual double eval(Eigen::Vector3d const &pos) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual double eval(Eigen::Vector3d const &pos) const = 0; }; /** * @brief Freestream potential - * @authors Adrien Crovato */ class FLOW_API F0PsPhiInf : public F0Ps { @@ -73,12 +47,9 @@ class FLOW_API F0PsPhiInf : public F0Ps public: F0PsPhiInf(int _nDim, double _alpha, double _beta = 0); - void update(double _alpha, double _beta = 0); + virtual void update(double _alpha) override; -#ifndef SWIG virtual double eval(Eigen::Vector3d const &pos) const override; - virtual void write(std::ostream &out) const override; -#endif }; } // namespace flow diff --git a/flow/src/wF1Ct.cpp b/flow/src/wF1Ct.cpp new file mode 100644 index 0000000000000000000000000000000000000000..57769dfe5787df52fb8bd92f7877c05843ebca5f --- /dev/null +++ b/flow/src/wF1Ct.cpp @@ -0,0 +1,154 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "wF1Ct.h" +using namespace tbox; +using namespace flow; + +F1CtDrag::F1CtDrag(int _nDim, double _sRef, double alpha, double _beta) : F1Ct(), nDim(_nDim), sRef(_sRef), beta(_beta) +{ + // Sanity check + if (nDim != 2 && nDim != 3) + { + std::stringstream err; + err << "F1CtDrag::nDim should be equal to 2 or 3 but " << nDim << " was given!\n"; + throw std::runtime_error(err.str()); + } + if (nDim != 3 && beta != 0) + throw std::runtime_error("F1CtDrag::beta can be nonzero only for 3D cases!\n"); + // Set the velocity + this->update(alpha); +} +/** + * @brief Update the direction and its gradient + */ +void F1CtDrag::update(double alpha) +{ + if (nDim == 2) + { + v = Eigen::Vector3d(cos(alpha), sin(alpha), 0) / sRef; + dv = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / sRef; + } + else + { + v = Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)) / sRef; + dv = Eigen::Vector3d(-sin(alpha) * cos(beta), sin(beta), cos(alpha) * cos(beta)) / sRef; + } +} +/** + * @brief Return the direction + */ +Eigen::Vector3d F1CtDrag::eval() const +{ + return v; +} +/** + * @brief Return the gradient of the direction with respect to angle of attack + */ +Eigen::Vector3d F1CtDrag::evalGrad() const +{ + return dv; +} + +F1CtSide::F1CtSide(int _nDim, double _sRef, double alpha, double _beta) : F1Ct(), nDim(_nDim), sRef(_sRef), beta(_beta) +{ + // Sanity check + if (nDim != 2 && nDim != 3) + { + std::stringstream err; + err << "F1CtSide::nDim should be equal to 2 or 3 but " << nDim << " was given!\n"; + throw std::runtime_error(err.str()); + } + if (nDim != 3 && beta != 0) + throw std::runtime_error("F1CtSide::beta can be nonzero only for 3D cases!\n"); + // Set the velocity + this->update(alpha); +} +/** + * @brief Update the direction and its gradient + */ +void F1CtSide::update(double alpha) +{ + if (nDim == 2) + { + v = Eigen::Vector3d(0, 0, 0); + dv = Eigen::Vector3d(0, 0, 0); + } + else + { + v = Eigen::Vector3d(-cos(alpha) * sin(beta), cos(beta), -sin(alpha) * sin(beta)) / sRef; + dv = Eigen::Vector3d(sin(alpha) * sin(beta), cos(beta), -cos(alpha) * sin(beta)) / sRef; + } +} +/** + * @brief Return the direction + */ +Eigen::Vector3d F1CtSide::eval() const +{ + return v; +} +/** + * @brief Return the gradient of the direction with respect to angle of attack + */ +Eigen::Vector3d F1CtSide::evalGrad() const +{ + return dv; +} + +F1CtLift::F1CtLift(int _nDim, double _sRef, double alpha, double _beta) : F1Ct(), nDim(_nDim), sRef(_sRef), beta(_beta) +{ + // Sanity check + if (nDim != 2 && nDim != 3) + { + std::stringstream err; + err << "F1CtLift::nDim should be equal to 2 or 3 but " << nDim << " was given!\n"; + throw std::runtime_error(err.str()); + } + if (nDim != 3 && beta != 0) + throw std::runtime_error("F1CtLift::beta can be nonzero only for 3D cases!\n"); + // Set the velocity + this->update(alpha); +} +/** + * @brief Update the direction and its gradient + */ +void F1CtLift::update(double alpha) +{ + if (nDim == 2) + { + v = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / sRef; + dv = Eigen::Vector3d(-cos(alpha), -sin(alpha), 0) / sRef; + } + else + { + v = Eigen::Vector3d(-sin(alpha), 0, cos(alpha)) / sRef; + dv = Eigen::Vector3d(-cos(alpha), 0, -sin(alpha)) / sRef; + } +} +/** + * @brief Return the direction + */ +Eigen::Vector3d F1CtLift::eval() const +{ + return v; +} +/** + * @brief Return the gradient of the direction with respect to angle of attack + */ +Eigen::Vector3d F1CtLift::evalGrad() const +{ + return dv; +} diff --git a/flow/src/wF1Ct.h b/flow/src/wF1Ct.h new file mode 100644 index 0000000000000000000000000000000000000000..7af9fe21b4496232f6cadde46ce277bdc90d2559 --- /dev/null +++ b/flow/src/wF1Ct.h @@ -0,0 +1,99 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef WF1CT_H +#define WF1CT_H + +#include "flow.h" +#include "wObserver.h" +#include <Eigen/Dense> + +namespace flow +{ + +/** + * @brief Constant vector function + */ +class FLOW_API F1Ct : public fwk::Observer +{ +public: + F1Ct() : Observer() {} + virtual ~F1Ct() {} + + virtual Eigen::Vector3d eval() const = 0; + virtual Eigen::Vector3d evalGrad() const = 0; +}; + +/** + * @brief Drag direction + */ +class FLOW_API F1CtDrag : public F1Ct +{ + int nDim; ///< number of dimensions + double sRef; ///< reference area + double beta; ///< angle of sideslip + Eigen::Vector3d v; ///< drag direction unit vector + Eigen::Vector3d dv; ///< gradient of drag direction unit vector (wrt alpha) + +public: + F1CtDrag(int _nDim, double _sRef, double alpha, double _beta = 0); + void update(double alpha); + + virtual Eigen::Vector3d eval() const override; + virtual Eigen::Vector3d evalGrad() const override; +}; + +/** + * @brief Sideforce direction + */ +class FLOW_API F1CtSide : public F1Ct +{ + int nDim; ///< number of dimensions + double sRef; ///< reference area + double beta; ///< angle of sideslip + Eigen::Vector3d v; ///< sideforce direction unit vector + Eigen::Vector3d dv; ///< gradient of sideforce direction unit vector (wrt alpha) + +public: + F1CtSide(int _nDim, double _sRef, double alpha, double _beta = 0); + void update(double alpha); + + virtual Eigen::Vector3d eval() const override; + virtual Eigen::Vector3d evalGrad() const override; +}; + +/** + * @brief Lift direction + */ +class FLOW_API F1CtLift : public F1Ct +{ + int nDim; ///< number of dimensions + double sRef; ///< reference area + double beta; ///< angle of sideslip + Eigen::Vector3d v; ///< lift direction unit vector + Eigen::Vector3d dv; ///< gradient of lift direction unit vector (wrt alpha) + +public: + F1CtLift(int _nDim, double _sRef, double alpha, double _beta = 0); + void update(double alpha); + + virtual Eigen::Vector3d eval() const override; + virtual Eigen::Vector3d evalGrad() const override; +}; + +} // namespace flow + +#endif //WF1CT_H diff --git a/flow/src/wF1El.cpp b/flow/src/wF1El.cpp index a65f35c86e7aedf64c7c30ff0c140394269c0d0c..d3d41fb781d4e8f472a09b7e70b092921631cde8 100644 --- a/flow/src/wF1El.cpp +++ b/flow/src/wF1El.cpp @@ -18,7 +18,7 @@ using namespace tbox; using namespace flow; -F1ElVi::F1ElVi(int _nDim, double alpha, double beta) : Fct1(), nDim(_nDim) +F1ElVi::F1ElVi(int _nDim, double alpha, double _beta) : F1El(), nDim(_nDim), beta(_beta) { // Sanity check if (nDim != 2 && nDim != 3) @@ -30,13 +30,12 @@ F1ElVi::F1ElVi(int _nDim, double alpha, double beta) : Fct1(), nDim(_nDim) if (nDim != 3 && beta != 0) throw std::runtime_error("F1ElVi::beta can be nonzero only for 3D cases!\n"); // Set the velocity - this->update(alpha, beta); + this->update(alpha); } /** * @brief Update the velocity vector - * @authors Adrien Crovato */ -void F1ElVi::update(double alpha, double beta) +void F1ElVi::update(double alpha) { if (nDim == 2) { @@ -51,7 +50,6 @@ void F1ElVi::update(double alpha, double beta) } /** * @brief Return the freestream velocity - * @authors Adrien Crovato */ Eigen::Vector3d F1ElVi::eval(Element const &e, std::vector<double> const &u, size_t k) const { @@ -60,11 +58,7 @@ Eigen::Vector3d F1ElVi::eval(Element const &e, std::vector<double> const &u, siz /** * @brief Return the gradient of freestream velocity with respect to angle of attack */ -Eigen::Vector3d F1ElVi::evalD(Element const &e, std::vector<double> const &u, size_t k) const +Eigen::Vector3d F1ElVi::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const { return dv; } -void F1ElVi::write(std::ostream &out) const -{ - out << "F1ElVi"; -} \ No newline at end of file diff --git a/flow/src/wF1El.h b/flow/src/wF1El.h index 7d86971217a83bf1705fa96dc3911fc1b48f3523..482b13e3d2b046fc7bc831205e07cbf3b6f5c866 100644 --- a/flow/src/wF1El.h +++ b/flow/src/wF1El.h @@ -18,36 +18,46 @@ #define WF1EL_H #include "flow.h" -#include "wFct1.h" +#include "wObserver.h" #include <vector> #include <Eigen/Dense> namespace flow { +/** + * @brief Vector function to be integrated over an element + */ +class FLOW_API F1El : public fwk::Observer +{ +public: + F1El() : Observer() {} + virtual ~F1El() {} + + virtual Eigen::Vector3d eval(tbox::Element const &e, std::vector<double> const &u, + size_t k) const = 0; + virtual Eigen::Vector3d evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const = 0; +}; + /** * @brief Freestream velocity - * @authors Adrien Crovato */ -class FLOW_API F1ElVi : public tbox::Fct1 +class FLOW_API F1ElVi : public F1El { int nDim; ///< number of dimensions + double beta; ///< angle of sideslip Eigen::Vector3d v; ///< freestream velocity vector Eigen::Vector3d dv; ///< freestream velocity vector gradient (wrt alpha) public: - F1ElVi(int _nDim, double alpha, double beta = 0); - void update(double alpha, double beta = 0); + F1ElVi(int _nDim, double alpha, double _beta = 0); + virtual void update(double alpha) override; -#ifndef SWIG - virtual Eigen::Vector3d eval(tbox::Element const &e, - std::vector<double> const &u, + virtual Eigen::Vector3d eval(tbox::Element const &e, std::vector<double> const &u, size_t k) const override; - virtual Eigen::Vector3d evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override; - virtual void write(std::ostream &out) const override; -#endif + virtual Eigen::Vector3d evalGrad(tbox::Element const &e, std::vector<double> const &u, + size_t k) const override; }; } // namespace flow diff --git a/flow/src/wFreestream.cpp b/flow/src/wFreestream.cpp index b319b5fe9fb61464db01269425bdde1d9056ef78..3ade264b298c192dd06236050e9af6dd57f987dd 100644 --- a/flow/src/wFreestream.cpp +++ b/flow/src/wFreestream.cpp @@ -15,17 +15,23 @@ */ #include "wFreestream.h" -#include "wProblem.h" #include "wTag.h" using namespace tbox; using namespace flow; -Freestream::Freestream(std::shared_ptr<MshData> _msh, int no, Fct1 &_f) : Group(_msh, no), f(_f) +Freestream::Freestream(std::shared_ptr<MshData> _msh, int no, int dim, double alpha, double beta) : Group(_msh, no) { + f = new F1ElVi(dim, alpha, beta); } -Freestream::Freestream(std::shared_ptr<MshData> _msh, std::string const &name, Fct1 &_f) : Group(_msh, name), f(_f) +Freestream::Freestream(std::shared_ptr<MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Group(_msh, name) { + f = new F1ElVi(dim, alpha, beta); +} +Freestream::~Freestream() +{ + delete f; + std::cout << "~Freestream()\n"; } void Freestream::write(std::ostream &out) const diff --git a/flow/src/wFreestream.h b/flow/src/wFreestream.h index 1725a2543da0bfc41b19b2c57b323d8bc471cd19..96a7871931d359a3f923054628dc1f4a6cc8e547 100644 --- a/flow/src/wFreestream.h +++ b/flow/src/wFreestream.h @@ -19,24 +19,24 @@ #include "flow.h" #include "wGroup.h" -#include "wFct1.h" +#include "wF1EL.h" namespace flow { /** - * @brief Manage freestream (Freestream) boundary condtion + * @brief Manage freestream (Neumann) boundary condtion * @authors Adrien Crovato */ class FLOW_API Freestream : public tbox::Group { public: #ifndef SWIG - tbox::Fct1 &f; ///< vector valued function to compute flux + F1ElVi *f; ///< vector valued function to compute flux #endif - Freestream(std::shared_ptr<tbox::MshData> _msh, int no, tbox::Fct1 &_f); - Freestream(std::shared_ptr<tbox::MshData> _msh, std::string const &name, tbox::Fct1 &_f); - virtual ~Freestream() { std::cout << "~Freestream()\n"; } + Freestream(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta = 0.0); + Freestream(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta = 0.0); + virtual ~Freestream(); #ifndef SWIG virtual void write(std::ostream &out) const override; diff --git a/flow/src/wFreestreamResidual.cpp b/flow/src/wFreestreamResidual.cpp index b2678b4501a4f57de19cc47146fbd58976537b30..04c8a24783714589cc0239406dde8320f11eeb74 100644 --- a/flow/src/wFreestreamResidual.cpp +++ b/flow/src/wFreestreamResidual.cpp @@ -36,7 +36,7 @@ Eigen::VectorXd FreestreamResidual::build(Element const &e, std::vector<double> Gauss &gauss = cache.getVGauss(); // Normal velocity - double vn = fs.f.eval(e, phi, 0).dot(e.getNormal()); + double vn = fs.f->eval(e, phi, 0).dot(e.getNormal()); // Integral of mass flux (density times normal velocity) through freestream boundary Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size()); @@ -57,7 +57,7 @@ Eigen::VectorXd FreestreamResidual::buildGradientAoa(tbox::Element const &e, std Gauss &gauss = cache.getVGauss(); // Gradient of normal velocity - double dvn = fs.f.evalD(e, phi, 0).dot(e.getNormal()); + double dvn = fs.f->evalGrad(e, phi, 0).dot(e.getNormal()); // Build mass flux gradient Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size()); diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp index 1113e3d8ed0702e7d9bae259ed983873fefdc6b1..401de65a83525ffb94201ebdaad9ad50b91fbf77 100644 --- a/flow/src/wLoadFunctional.cpp +++ b/flow/src/wLoadFunctional.cpp @@ -15,7 +15,7 @@ */ #include "wLoadFunctional.h" -#include "wFct0.h" +#include "wF0El.h" #include "wElement.h" #include "wCache.h" @@ -29,7 +29,7 @@ using namespace flow; * * b = \int psi * cp * n dS */ -Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp) +Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::vector<double> const &phi, F0El const &cp) { // Get pre-computed values Cache &cache = e.getVCache(); @@ -55,14 +55,14 @@ Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std:: * * A = \int psi * dcp * n dS */ -Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp) +Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element const &eV, std::vector<double> const &phi, F0El const &cp) { // Get pre-computed values Cache &cache = e.getVCache(); Gauss &gauss = cache.getVGauss(); // Gradient of pressure force coefficient in associated volume - Eigen::MatrixXd dfcp = -e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0); // "-" because force on body not on fluid + Eigen::MatrixXd dfcp = -e.getNormal() * cp.evalGrad(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0); // "-" because force on body not on fluid // Integrate cp on the surface Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size()); @@ -84,7 +84,7 @@ Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element cons * + \int psi * cp * dn * dS * + \int psi * cp * n * ddS */ -std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp, int nDim) +std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(Element const &e, Element const &eV, std::vector<double> const &phi, F0El const &cp, int nDim) { // Get pre-computed values Cache &cache = e.getVCache(); @@ -94,7 +94,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E double cP = -cp.eval(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid) Eigen::Vector3d nrm = e.getNormal(); // unit normal // Gradients - double dCp = -cp.evalD(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid) + double dCp = -cp.evalGrad(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid) Eigen::VectorXd dPhi = eV.computeGradient(phi, 0); // potential Eigen::MatrixXd const &iJ = eV.getJinv(0); // inverse Jacobian std::vector<Eigen::MatrixXd> const &dJ = eV.getGradJ(0); // Jacobian diff --git a/flow/src/wLoadFunctional.h b/flow/src/wLoadFunctional.h index 4a1a25a8ca3497dc18c1210d49b61dd20cb4c425..105fe2b7fcac05800630c59ce8fb54b7b38c43a9 100644 --- a/flow/src/wLoadFunctional.h +++ b/flow/src/wLoadFunctional.h @@ -32,10 +32,10 @@ class FLOW_API LoadFunctional { public: // Solver - static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, tbox::Fct0 const &cp); + static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, F0El const &cp); // Adjoint - static Eigen::MatrixXd buildGradientFlow(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, tbox::Fct0 const &cp); - static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, tbox::Fct0 const &cp, int nDim); + static Eigen::MatrixXd buildGradientFlow(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, F0El const &cp); + static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, F0El const &cp, int nDim); }; } // namespace flow diff --git a/flow/src/wMedium.cpp b/flow/src/wMedium.cpp index 9d4fecc2c3dda3016409ee22f5992aaaed0d80c5..c03d9a81d65cbde29ac3a2d00067a4b5c8397c37 100644 --- a/flow/src/wMedium.cpp +++ b/flow/src/wMedium.cpp @@ -21,27 +21,54 @@ using namespace tbox; using namespace flow; -Medium::Medium(std::shared_ptr<MshData> _msh, int no, Fct0 &_rho, Fct0 &_mach, - Fct0 &_cP, F0PsPhiInf &_phiInf) : Group(_msh, no), rho(_rho), - mach(_mach), cP(_cP), phiInf(_phiInf) +Medium::Medium(std::shared_ptr<MshData> _msh, int no, + double mInf, int dim, double alpha, double beta) : Group(_msh, no) { + initFun(mInf, dim, alpha, beta); createMap(); std::cout << "Fluid is " << *tag << "with " << nCnt << " nodes." << std::endl << std::endl; } Medium::Medium(std::shared_ptr<MshData> _msh, std::string const &name, - Fct0 &_rho, Fct0 &_mach, Fct0 &_cP, F0PsPhiInf &_phiInf) : Group(_msh, name), rho(_rho), mach(_mach), cP(_cP), phiInf(_phiInf) + double mInf, int dim, double alpha, double beta) : Group(_msh, name) { + initFun(mInf, dim, alpha, beta); createMap(); std::cout << "Fluid is " << *tag << "with " << nCnt << " nodes." << std::endl << std::endl; } +Medium::~Medium() +{ + delete rho; + delete mach; + delete cP; + std::cout << "~Medium()\n"; +} + +/** + * Initialize functions + */ +void Medium::initFun(double mInf, int dim, double alpha, double beta) +{ + if (mInf == 0) + { + rho = new F0ElRhoL(); + mach = new F0ElMachL(); + cP = new F0ElCpL(); + } + else + { + rho = new F0ElRho(mInf); + mach = new F0ElMach(mInf); + cP = new F0ElCp(mInf); + } + phiInf = new F0PsPhiInf(dim, alpha, beta); +} /** * @brief Create data structure - * @authors Adrien Crovato */ void Medium::createMap() { diff --git a/flow/src/wMedium.h b/flow/src/wMedium.h index f84b815d56be1c01399978386196adc8f65c558a..7bfbe69116b77838639052e95e6464ecce47b176 100644 --- a/flow/src/wMedium.h +++ b/flow/src/wMedium.h @@ -19,8 +19,9 @@ #include "flow.h" #include "wGroup.h" -#include "wFct0.h" +#include "wF0El.h" #include "wF0Ps.h" +#include <map> namespace flow { @@ -31,24 +32,26 @@ namespace flow */ class FLOW_API Medium : public tbox::Group { +private: + size_t nCnt; ///< number of nodes public: std::map<tbox::Node *, std::vector<tbox::Element *>> neMap; ///< map each elements connected to a node - size_t nCnt; ///< number of nodes #ifndef SWIG - tbox::Fct0 ρ ///< density - tbox::Fct0 &mach; ///< Mach number - tbox::Fct0 &cP; ///< pressure coefficient - flow::F0PsPhiInf &phiInf; ///< freestream potential + F0El *rho; ///< density + F0El *mach; ///< Mach number + F0El *cP; ///< pressure coefficient + F0PsPhiInf *phiInf; ///< freestream potential #endif std::vector<std::pair<tbox::Element *, std::vector<tbox::Element *>>> adjMap; ///< map of adjacent elements - Medium(std::shared_ptr<tbox::MshData> _msh, int no, tbox::Fct0 &_rho, tbox::Fct0 &_mach, tbox::Fct0 &_cP, flow::F0PsPhiInf &_phiInf); - Medium(std::shared_ptr<tbox::MshData> _msh, std::string const &name, tbox::Fct0 &_rho, tbox::Fct0 &_mach, tbox::Fct0 &_cP, flow::F0PsPhiInf &_phiInf); - virtual ~Medium() { std::cout << "~Medium()\n"; } + Medium(std::shared_ptr<tbox::MshData> _msh, int no, double mInf, int dim, double alpha, double beta); + Medium(std::shared_ptr<tbox::MshData> _msh, std::string const &name, double mInf, int dim, double alpha, double beta); + virtual ~Medium(); #ifndef SWIG + void initFun(double mInf, int dim, double alpha, double beta); void createMap(); virtual void write(std::ostream &out) const override; #endif diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp index ec0a579a6825e21c5ae398317124b26d45443904..4bff90087c57ae4ba84f15ca1e9a8d2b196a9e04 100644 --- a/flow/src/wNewton.cpp +++ b/flow/src/wNewton.cpp @@ -287,7 +287,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J) } } // Assembly (supersonic) - double mach = fluid->mach.eval(*e, phi, 0); + double mach = fluid->mach->eval(*e, phi, 0); if (mach > mCOv) { for (size_t ii = 0; ii < e->nodes.size(); ++ii) @@ -309,11 +309,10 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J) tms["1-Jkutta"].start(); for (auto wake : pbl->wBCs) { - Eigen::VectorXd vi = pbl->computeVi(); tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) { // Build Ku, Kl Eigen::MatrixXd Aue, Ale; - std::tie(Aue, Ale) = WakeResidual::buildGradientFlow(*we, phi, vi); + std::tie(Aue, Ale) = WakeResidual::buildGradientFlow(*we, phi); // Assembly tbb::spin_mutex::scoped_lock lock(mutex); for (size_t ii = 0; ii < we->nRow; ++ii) @@ -431,7 +430,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R) // Apply blowing (transpiration) BCs for (auto bBC : pbl->bBCs) { - tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, Fct0C *> p) { + tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, F0ElC *> p) { // Build elementary flux vector Eigen::VectorXd be = BlowingResidual::build(*p.first, phi, *p.second); @@ -451,11 +450,10 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R) // Nishida Thesis, pp. 29-30, Eq. 2.18 for (auto wake : pbl->wBCs) { - Eigen::VectorXd vi = pbl->computeVi(); tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) { // Build bu, bl Eigen::VectorXd bue, ble; - std::tie(bue, ble) = WakeResidual::build(*we, phi, vi); + std::tie(bue, ble) = WakeResidual::build(*we, phi); // Assembly tbb::spin_mutex::scoped_lock lock(mutex); for (size_t ii = 0; ii < we->nRow; ++ii) diff --git a/flow/src/wPicard.cpp b/flow/src/wPicard.cpp index 8b66292b19a49e4be04a2355e71ddbb91718ea9d..4a3a890b5a3faacc877f8b6de6d1e7b906b40c3a 100644 --- a/flow/src/wPicard.cpp +++ b/flow/src/wPicard.cpp @@ -226,11 +226,10 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector< tms["1-Akutta"].start(); for (auto wake : pbl->wBCs) { - Eigen::VectorXd vi = pbl->computeVi(); tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) { // Build Ku, Kl Eigen::MatrixXd Aue, Ale; - std::tie(Aue, Ale) = WakeResidual::buildFixed(*we, phi, vi); + std::tie(Aue, Ale) = WakeResidual::buildFixed(*we, phi); // Assembly tbb::spin_mutex::scoped_lock lock(mutex); for (size_t ii = 0; ii < we->nRow; ++ii) @@ -293,7 +292,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector< // Apply blowing (transpiration) BCs for (auto bBC : pbl->bBCs) { - tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, Fct0C *> p) { + tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, F0ElC *> p) { // Build elementary flux vector Eigen::VectorXd be = BlowingResidual::build(*p.first, phi, *p.second); diff --git a/flow/src/wPotentialResidual.cpp b/flow/src/wPotentialResidual.cpp index 1e683747272a8b1da7f599e23c1e482dd2d7f0e4..9122b663326c241625b573d73a0e29c8140ecc17 100644 --- a/flow/src/wPotentialResidual.cpp +++ b/flow/src/wPotentialResidual.cpp @@ -42,7 +42,7 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub // Shape functions gradient Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k); // Elementary stiffness matrix - A += dSf.transpose() * dSf * fluid.rho.eval(e, phi, k) * gauss.getW(k) * e.getDetJ(k); + A += dSf.transpose() * dSf * fluid.rho->eval(e, phi, k) * gauss.getW(k) * e.getDetJ(k); } return A; } @@ -51,7 +51,6 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub * @brief Build the residual vector, on one element * * b = \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi dV - * @todo remove muC, mCO + mu as function? */ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO) { @@ -62,14 +61,14 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st // Subsonic contribution Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size()); for (size_t k = 0; k < gauss.getN(); ++k) - b += fluid.rho.eval(e, phi, k) * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k); + 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 mach = fluid.mach->eval(e, phi, 0); if (mach > mCO) { double mu = muC * (1 - (mCO * mCO) / (mach * mach)); // switching function - double rhoU = fluid.rho.eval(eU, phi, 0); // density on upwind element + double rhoU = fluid.rho->eval(eU, phi, 0); // density on upwind element // scale 1st term b *= 1 - mu; // contribution of upwind element @@ -87,7 +86,6 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st * + \int mu*drhoU * grad_phi * grad_psi dV * + \int ( (1-mu)*rho + mu*rhoU ) * dgrad_phi * grad_psi dV * + \int (rhoU - rho)*dmu * grad_phi * grad_psi dV - * @todo remove muC, mCO + mu as function? */ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO) { @@ -106,13 +104,13 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo Eigen::VectorXd dPhi = e.computeGradient(phi, k); // rho * grad_phi*grad_psi + drho * grad_phi*grad_psi - A1 += fluid.rho.eval(e, phi, k) * dSf.transpose() * dSf * wdj; - A1 += fluid.rho.evalD(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj; + A1 += fluid.rho->eval(e, phi, k) * dSf.transpose() * dSf * wdj; + A1 += fluid.rho->evalGrad(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj; } // Supersonic contribution Eigen::MatrixXd A2; - double mach = fluid.mach.eval(e, phi, 0); + double mach = fluid.mach->eval(e, phi, 0); if (mach > mCO) { // switching function and gradient @@ -120,8 +118,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo double dmu = 2 * muC * mCO * mCO / (mach * mach * mach); // density and gradients on upwind elements Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0); - double rhoU = fluid.rho.eval(eU, phi, 0); - double dRhoU = fluid.rho.evalD(eU, phi, 0); + double rhoU = fluid.rho->eval(eU, phi, 0); + double dRhoU = fluid.rho->evalGrad(eU, phi, 0); Eigen::MatrixXd const &dSfU = eU.getJinv(0) * eU.getVCache().getDsf(0); // upwinding @@ -139,7 +137,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo A2 += mu * dRhoU * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj; // mu*rhoU*grad_phi*grad_psi + dmu*(rhoU-rho)*grad_phi*grad_psi A1 += mu * rhoU * dSf.transpose() * dSf * wdj; - A1 += dmu * (rhoU - fluid.rho.eval(e, phi, k)) * fluid.mach.evalD(e, phi, k) * dSf.transpose() * dPhi * dPhi.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; } } return std::make_tuple(A1, A2); @@ -155,7 +153,6 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo * + \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * dgrad_psi dV * + \int (rhoU - rho)*dmu * grad_phi * grad_psi dV * + \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi ddV - * @todo remove muC, mCO + mu as function? */ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, int nDim, double muC, double mCO) { @@ -169,8 +166,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes { // Density and gradients Eigen::VectorXd dPhi = e.computeGradient(phi, k); - double rho = fluid.rho.eval(e, phi, k); - double dRho = fluid.rho.evalD(e, phi, k); + double rho = fluid.rho->eval(e, phi, k); + double dRho = fluid.rho->evalGrad(e, phi, k); Eigen::MatrixXd const &dSf = cache.getDsf(k); Eigen::MatrixXd const &iJ = e.getJinv(k); // Jacobian gradients @@ -193,7 +190,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes } // Supersonic contribution Eigen::MatrixXd A2; - double mach = fluid.mach.eval(e, phi, 0); + double mach = fluid.mach->eval(e, phi, 0); if (mach > mCO) { // switching function and gradient @@ -201,8 +198,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes double dmu = 2 * muC * mCO * mCO / (mach * mach * mach); // upwind density and gradient Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0); - double rhoU = fluid.rho.eval(eU, phi, 0); - double dRhoU = fluid.rho.evalD(eU, phi, 0); + double rhoU = fluid.rho->eval(eU, phi, 0); + double dRhoU = fluid.rho->evalGrad(eU, phi, 0); Eigen::MatrixXd const &iJU = eU.getJinv(0); std::vector<Eigen::MatrixXd> const &dJU = eU.getGradJ(0); @@ -213,8 +210,8 @@ 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.evalD(e, 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 diff --git a/flow/src/wPotentialResidual.h b/flow/src/wPotentialResidual.h index 62dde67b704aab847e97578315053d03b97c432e..7972271854c352ceaea3b01679ff816d35e91b72 100644 --- a/flow/src/wPotentialResidual.h +++ b/flow/src/wPotentialResidual.h @@ -27,6 +27,7 @@ namespace flow /** * @brief Formulation of nonlinear potential equation residuals + * @todo mu as F0ElMu? * @todo split stabilization terms from subsonic terms? */ class FLOW_API PotentialResidual diff --git a/flow/src/wProblem.cpp b/flow/src/wProblem.cpp index 7d3d069fe05aa57a3b2cc416cff3b621e52026de..46bc8ad6b357156b38a1b3e32193d0788d3220fe 100644 --- a/flow/src/wProblem.cpp +++ b/flow/src/wProblem.cpp @@ -29,96 +29,29 @@ using namespace flow; Problem::Problem(std::shared_ptr<MshData> _msh, int dim, double aoa, double aos, double minf, double sref, double cref, double xref, double yref, double zref) : msh(_msh), nDim(dim), alpha(aoa), - beta(aos), M_inf(minf), S_ref(sref), - c_ref(cref) + beta(aos), M_inf(minf), + S_ref(sref), c_ref(cref), + x_ref(xref, yref, zref), + dirD(dim, sref, aoa, aos), + dirS(dim, sref, aoa, aos), + dirL(dim, sref, aoa, aos) { - x_ref(0) = xref; - x_ref(1) = yref; - x_ref(2) = zref; -} - -/** - * @brief Compute drag, sideforce and lift normalized directions - * @todo use F1ElVi - */ -std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> Problem::computeDirections() -{ - Eigen::Vector3d dirD, dirS, dirL; - if (nDim == 2) - { - dirD = Eigen::Vector3d(cos(alpha), sin(alpha), 0) / S_ref; - dirS = Eigen::Vector3d(0, 0, 0); - dirL = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / S_ref; - } - else - { - dirD = Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)) / S_ref; - dirS = Eigen::Vector3d(-cos(alpha) * sin(beta), cos(beta), -sin(alpha) * sin(beta)) / S_ref; - dirL = Eigen::Vector3d(-sin(alpha), 0, cos(alpha)) / S_ref; - } - return std::make_tuple(dirD, dirS, dirL); -} - -/** - * @brief Compute gradient of drag, sideforce and lift normalized directions - * @todo use F1ElVi - * @todo only wrt alpha - */ -std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> Problem::computeGradientDirections() -{ - Eigen::Vector3d dirD, dirS, dirL; - if (nDim == 2) - { - dirD = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / S_ref; - dirS = Eigen::Vector3d(0, 0, 0); - dirL = Eigen::Vector3d(-cos(alpha), -sin(alpha), 0) / S_ref; - } - else - { - dirD = Eigen::Vector3d(-sin(alpha) * cos(beta), sin(beta), cos(alpha) * cos(beta)) / S_ref; - dirS = Eigen::Vector3d(sin(alpha) * sin(beta), cos(beta), -cos(alpha) * sin(beta)) / S_ref; - dirL = Eigen::Vector3d(-cos(alpha), 0, -sin(alpha)) / S_ref; - } - return std::make_tuple(dirD, dirS, dirL); -} - -/** - * @brief Compute freestream velocity - * @todo use F1ElVi - */ -Eigen::VectorXd Problem::computeVi() -{ - if (nDim == 2) - return Eigen::Vector2d(cos(alpha), sin(alpha)); - else - return Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)); -} - -/** - * @brief Compute freestream velocity - * @todo use F1ElVi - * @todo only wrt alpha - */ -Eigen::VectorXd Problem::computeGradientVi() -{ - if (nDim == 2) - return Eigen::Vector2d(-sin(alpha), cos(alpha)); - else - return Eigen::Vector3d(-sin(alpha) * cos(beta), 0, cos(alpha) * cos(beta)); + obs.push_back(&dirD); + obs.push_back(&dirS); + obs.push_back(&dirL); } /** * @brief Set the fluid medium - * @authors Adrien Crovato */ void Problem::set(std::shared_ptr<Medium> m) { medium = m; + obs.push_back(m->phiInf); } /** * @brief Add a boundary surface - * @authors Adrien Crovato */ void Problem::add(std::shared_ptr<Boundary> b) { @@ -127,34 +60,33 @@ void Problem::add(std::shared_ptr<Boundary> b) /** * @brief Add initial condition - * @authors Adrien Crovato */ void Problem::set(std::shared_ptr<Initial> i) { iIC = i; + obs.push_back(i->f); } /** * @brief Add Dirichlet boundary condition - * @authors Adrien Crovato */ void Problem::add(std::shared_ptr<Dirichlet> d) { dBCs.push_back(d); + obs.push_back(d->f); } /** * @brief Add freestream (Neumann) boundary condition - * @authors Adrien Crovato */ void Problem::add(std::shared_ptr<Freestream> f) { fBCs.push_back(f); + obs.push_back(f->f); } /** * @brief Add wake boundary condition - * @authors Adrien Crovato */ void Problem::add(std::shared_ptr<Wake> w) { @@ -163,7 +95,6 @@ void Problem::add(std::shared_ptr<Wake> w) /** * @brief Add Kutta condition - * @authors Adrien Crovato */ void Problem::add(std::shared_ptr<Kutta> k) { @@ -172,13 +103,22 @@ void Problem::add(std::shared_ptr<Kutta> k) /** * @brief Add Blowing boundary condition - * @authors Adrien Crovato */ void Problem::add(std::shared_ptr<Blowing> b) { bBCs.push_back(b); } +/** + * @brief Update the angle of attack + */ +void Problem::update(double aoa) +{ + alpha = aoa; + for (auto o : obs) + o->update(aoa); +} + /** * @brief Initialize the elements precomputed values */ @@ -186,7 +126,7 @@ void Problem::initElems() { // Update volume for (auto e : medium->tag->elems) - e->initValues(true); + e->initValues(true, 1); // Update surface (boundaries) for (auto surf : bnds) for (auto e : surf->groups[0]->tag->elems) @@ -212,7 +152,6 @@ void Problem::initElems() /** * @brief Initialize the elements precomputed gradients - * @authors Adrien Crovato */ void Problem::initGradElems() { @@ -244,12 +183,11 @@ void Problem::initGradElems() void Problem::pin() { if (dBCs.empty()) - this->add(std::make_shared<Dirichlet>(this->msh, this->fBCs.front()->tag->no, this->medium->phiInf, true)); + this->add(std::make_shared<Dirichlet>(msh, fBCs.front()->tag->no, nDim, alpha, beta, true)); } /** * @brief Check that Problem is not empty and that elements are supported - * @authors Adrien Crovato */ void Problem::check() const { diff --git a/flow/src/wProblem.h b/flow/src/wProblem.h index 152aa1f267cc0530190b5c386f8ecf2b591ebac6..53f2e5b3024dfa7cb30e68b95f5b839f76a8640b 100644 --- a/flow/src/wProblem.h +++ b/flow/src/wProblem.h @@ -19,6 +19,7 @@ #include "flow.h" #include "wObject.h" +#include "wF1Ct.h" #include <memory> #include <iostream> #include <vector> @@ -30,42 +31,45 @@ namespace flow /** * @brief Manage the problem * - * Contains mesh, medium, boundaries and reference value + * Contains freestream definition, reference values and physical groups * @authors Adrien Crovato */ class FLOW_API Problem : public fwk::wSharedObject { public: std::shared_ptr<tbox::MshData> msh; ///< Mesh structure + int nDim; ///< Problem dimension + // Freestream + double alpha; ///< Angle of attack + double beta; ///< Angle of sideslip + double M_inf; ///< Mach number + // Reference values + double S_ref; ///< Reference surface + double c_ref; ///< Reference chord + Eigen::Vector3d x_ref; ///< Reference center point (for moment computation) #ifndef SWIG - std::shared_ptr<Medium> medium; ///< Fluid - std::vector<std::shared_ptr<Boundary>> bnds; ///< Boundaries - + // Loads direction + F1CtDrag dirD; // drag + F1CtSide dirS; // sideforce + F1CtLift dirL; // lift + // Physical groups + std::shared_ptr<Medium> medium; ///< Fluid + std::vector<std::shared_ptr<Boundary>> bnds; ///< Boundaries std::shared_ptr<Initial> iIC; ///< Initial condition std::vector<std::shared_ptr<Dirichlet>> dBCs; ///< Dirichlet boundary conditions std::vector<std::shared_ptr<Freestream>> fBCs; ///< Freestream boundary conditions std::vector<std::shared_ptr<Wake>> wBCs; ///< Wake boundary condition std::vector<std::shared_ptr<Kutta>> kSCs; ///< Kutta condition std::vector<std::shared_ptr<Blowing>> bBCs; ///< Blowing (transpiration) boundary condition + // Observers + std::vector<fwk::Observer *> obs; ///< classes depending on problem variables #endif - // Reference values - int nDim; ///< Problem dimension - double alpha; ///< Angle of attack - double beta; ///< Angle of sideslip - double M_inf; ///< Freestream Mach number - double S_ref; ///< Reference surface - double c_ref; ///< Reference chord - Eigen::Vector3d x_ref; ///< Reference center point (for moment computation) Problem(std::shared_ptr<tbox::MshData> _msh, int dim, double aoa, double aos, double minf, double sref, double cref, double xref, double yref, double zref); virtual ~Problem() { std::cout << "~Problem()\n"; } - std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> computeDirections(); - std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> computeGradientDirections(); - Eigen::VectorXd computeVi(); - Eigen::VectorXd computeGradientVi(); void set(std::shared_ptr<Medium> m); void add(std::shared_ptr<Boundary> b); void set(std::shared_ptr<Initial> i); @@ -74,6 +78,7 @@ public: void add(std::shared_ptr<Wake> w); void add(std::shared_ptr<Kutta> k); void add(std::shared_ptr<Blowing> b); + void update(double aoa); #ifndef SWIG void check() const; diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp index 7fae0493b0a2c12bdee06baf62efc0fe6f609eb6..991721ff36e3fb6d317562733df113c2997b054c 100644 --- a/flow/src/wSolver.cpp +++ b/flow/src/wSolver.cpp @@ -160,17 +160,11 @@ void Solver::save(MshExport *mshWriter, int n) */ void Solver::computeLoad() { - // Reset coefficients + // Reset coefficients and compute loads Cl = 0; Cd = 0; Cs = 0; Cm = 0; - - // Lift, drag and sideforce normalized directions - Eigen::Vector3d dirL, dirD, dirS; - std::tie(dirD, dirS, dirL) = pbl->computeDirections(); - - // Compute loads for (auto sur : pbl->bnds) { // Reset @@ -180,7 +174,7 @@ void Solver::computeLoad() for (auto e : sur->groups[0]->tag->elems) { // Build nodal pressure load - Eigen::VectorXd be = LoadFunctional::build(*e, *sur->svMap.at(e), phi, pbl->medium->cP); + Eigen::VectorXd be = LoadFunctional::build(*e, *sur->svMap.at(e), phi, *pbl->medium->cP); // Assembly for (size_t i = 0; i < e->nodes.size(); ++i) @@ -201,9 +195,9 @@ void Solver::computeLoad() sur->Cm += (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up } // rotate to flow direction - sur->Cd = Cf.dot(dirD); - sur->Cs = Cf.dot(dirS); - sur->Cl = Cf.dot(dirL); + sur->Cd = Cf.dot(pbl->dirD.eval()); + sur->Cs = Cf.dot(pbl->dirS.eval()); + sur->Cl = Cf.dot(pbl->dirL.eval()); // compute total Cl += sur->Cl; Cd += sur->Cd; @@ -223,7 +217,7 @@ void Solver::computeFlow() tbb::parallel_do(fluid->neMap.begin(), fluid->neMap.end(), [&](std::pair<Node *, std::vector<Element *>> p) { Node *n = p.first; // Perturbation potential - vPhi[n->row] = phi[n->row] - fluid->phiInf.eval(n->pos); + vPhi[n->row] = phi[n->row] - fluid->phiInf->eval(n->pos); // Velocity (averaged at nodes) U[n->row] = Eigen::Vector3d::Zero(); for (auto e : p.second) @@ -236,17 +230,17 @@ void Solver::computeFlow() // Density (averaged at nodes) rho[n->row] = 0.; for (auto e : p.second) - rho[n->row] += fluid->rho.eval(*e, phi, 0); + rho[n->row] += fluid->rho->eval(*e, phi, 0); rho[n->row] /= p.second.size(); // Mach number (averaged at nodes) M[n->row] = 0.; for (auto e : p.second) - M[n->row] += fluid->mach.eval(*e, phi, 0); + M[n->row] += fluid->mach->eval(*e, phi, 0); M[n->row] /= p.second.size(); // Pressure coefficient (averaged at nodes) Cp[n->row] = 0.; for (auto e : p.second) - Cp[n->row] += fluid->cP.eval(*e, phi, 0); + Cp[n->row] += fluid->cP->eval(*e, phi, 0); Cp[n->row] /= p.second.size(); }); diff --git a/flow/src/wWakeResidual.cpp b/flow/src/wWakeResidual.cpp index 12562c7ca3b69befdb65c2c0b6e4475e8ededd96..fcfd06b7f44f3748fb55c41c73433f6ebfe1e72d 100644 --- a/flow/src/wWakeResidual.cpp +++ b/flow/src/wWakeResidual.cpp @@ -30,8 +30,9 @@ using namespace flow; * A_u = \int (psi + grad_psi*v_inf) * v_u*grad_phi_u dS * A_l = \int (psi + grad_psi*v_inf) * v_l*grad_phi_l dS * NB: A_u and A_l will be assembled on same row DOF, but on different columns + * NB: v_inf = [1, 0, 0] */ -std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi) +std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElement const &we, std::vector<double> const &phi) { // Get shape functions and Gauss points Cache &sCache = we.surUpE->getVCache(); @@ -51,7 +52,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen for (size_t j = 0; j < nDim; ++j) dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first)); Eigen::VectorXd dSfp(nRow); - dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * vUpJ.transpose() * vi; + dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * vUpJ.transpose() * Eigen::Vector3d(1, 0, 0).block(0, 0, nDim, 1); // Velocity Eigen::RowVectorXd Aup = we.volUpE->computeGradient(phi, 0).transpose() * vUpJ * vUpDsf; Eigen::RowVectorXd Alw = we.volLwE->computeGradient(phi, 0).transpose() * vLwJ * vLwDsf; @@ -78,8 +79,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen * * b_u = \int (psi + grad_psi*v_inf) * (grad_phi_u)^2 dS * b_l = \int (psi + grad_psi*v_inf) * (grad_phi_l)^2 dS + * NB: v_inf = [1, 0, 0] */ -std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi) +std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement const &we, std::vector<double> const &phi) { // Get shape functions and Gauss points Cache &sCache = we.surUpE->getVCache(); @@ -95,7 +97,7 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con for (size_t j = 0; j < nDim; ++j) dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first)); Eigen::VectorXd dSfp(nRow); - dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * vi; + dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * Eigen::Vector3d(1, 0, 0).block(0, 0, nDim, 1); // Velocity squared Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0); Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0); @@ -124,12 +126,13 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con * * A_u = \int (psi + grad_psi*v_inf) * 2*(grad_phi_u)*dgrad_phi_u dS * A_l = \int (psi + grad_psi*v_inf) * 2*(grad_phi_l)*dgrad_phi_l dS + * NB: v_inf = [1, 0, 0] */ -std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientFlow(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi) +std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientFlow(WakeElement const &we, std::vector<double> const &phi) { // Build Eigen::MatrixXd Au, Al; - std::tie(Au, Al) = WakeResidual::buildFixed(we, phi, vi); + std::tie(Au, Al) = WakeResidual::buildFixed(we, phi); return std::make_tuple(2 * Au, 2 * Al); } @@ -140,8 +143,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientFlow(Wak * = \int dgrad_psi*v_inf * (grad_phi)^2 dS * + \int (psi + grad_psi*v_inf) * d(grad_phi)^2 dS * + \int (psi + grad_psi*v_inf) * (grad_phi)^2 ddS + * NB: v_inf = [1, 0, 0] */ -std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientMesh(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi) +std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientMesh(WakeElement const &we, std::vector<double> const &phi) { // Get shape functions and Gauss points Cache &sCache = we.surUpE->getVCache(); @@ -162,6 +166,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil size_t nSur = we.surUpE->nodes.size(); // Stabilization term + Eigen::VectorXd vi = Eigen::Vector3d(1, 0, 0).block(0, 0, nDim, 1); double sqSurf = sqrt(we.surUpE->getVol()); Eigen::MatrixXd dSf(nRow, nDim); for (size_t i = 0; i < nRow; ++i) diff --git a/flow/src/wWakeResidual.h b/flow/src/wWakeResidual.h index fc536849ef6c9931c1494d9ba7b937dfbe057c26..5bad9298cfcaf6d3d6fef052c53eae12015ab8b9 100644 --- a/flow/src/wWakeResidual.h +++ b/flow/src/wWakeResidual.h @@ -32,12 +32,12 @@ class FLOW_API WakeResidual { public: // Picard - static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildFixed(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi); + static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildFixed(WakeElement const &we, std::vector<double> const &phi); // Newton - static std::tuple<Eigen::VectorXd, Eigen::VectorXd> build(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi); - static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi); + static std::tuple<Eigen::VectorXd, Eigen::VectorXd> build(WakeElement const &we, std::vector<double> const &phi); + static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(WakeElement const &we, std::vector<double> const &phi); // Adjoint - static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi); + static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(WakeElement const &we, std::vector<double> const &phi); }; } // namespace flow diff --git a/flow/tests/adjoint.py b/flow/tests/adjoint.py index e0134b0ffb8b19a77454cc9fdec4c85f9c05b6b2..990fb526cc3ccc39cba4a8134250b8b2bc28e3f9 100644 --- a/flow/tests/adjoint.py +++ b/flow/tests/adjoint.py @@ -43,7 +43,6 @@ def main(): # define flow variables alpha = 2*np.pi/180 M_inf = 0.7 - M_crit = 5 # Squared critical Mach number (above which density is modified) c_ref = 1 dim = 2 @@ -58,7 +57,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') tms['pre'].stop() # solve problem diff --git a/flow/tests/adjoint3.py b/flow/tests/adjoint3.py index 670ed6bf2f404cfe70cd18a1687fa0173038fe71..c52cf8d7a9175e6060e9da34d6946cd54f35829a 100644 --- a/flow/tests/adjoint3.py +++ b/flow/tests/adjoint3.py @@ -43,7 +43,6 @@ def main(): # define flow variables alpha = 3*np.pi/180 M_inf = 0.3 - M_crit = 5 # Squared critical Mach number (above which density is modified) c_ref = 1 span = 1 dim = 3 @@ -59,7 +58,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') tms['pre'].stop() # solve problem diff --git a/flow/tests/bli.py b/flow/tests/bli.py index 73728691de5e6f4a4389359f8d77d57d3781439e..3c7b3584852f6b5a98062430f9b7841771b27615 100644 --- a/flow/tests/bli.py +++ b/flow/tests/bli.py @@ -55,7 +55,6 @@ def main(): alpha = 5*math.pi/180 U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1 M_inf = 0.5 - M_crit = 5 # Squared critical Mach number (above which density is modified) c_ref = 1 dim = 2 tol = 1e-4 # tolerance for the VII (usually one drag count) @@ -69,7 +68,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) + pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) tms['pre'].stop() # solve viscous problem diff --git a/flow/tests/cylinder.py b/flow/tests/cylinder.py index 2a0b51c433df7a11d08b6b69297445d11d0bea4a..1cd3a0d9d2f6686bf33549545aabee1640101795 100644 --- a/flow/tests/cylinder.py +++ b/flow/tests/cylinder.py @@ -40,7 +40,6 @@ def main(): alpha = 3*math.pi/180 U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1 M_inf = 0.1 - M_crit = 5 # Squared critical Mach number (above which density is modified) dim = 2 # define dimension and mesh size @@ -60,7 +59,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, l_ref, l_ref, 0., 0., 0., 'cylinder', te = 'te', dbc=True) + pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, l_ref, l_ref, 0., 0., 0., 'cylinder', te = 'te', dbc=True) tms['pre'].stop() # solve problem diff --git a/flow/tests/cylinder2D5.py b/flow/tests/cylinder2D5.py index 14f943dfd35bb2a56079f9bcbe7ab3ddbbfb9393..3999b7e4333422d842dc146e5c32ae11e674754a 100644 --- a/flow/tests/cylinder2D5.py +++ b/flow/tests/cylinder2D5.py @@ -40,7 +40,6 @@ def main(): alpha = 3*math.pi/180 U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1 M_inf = 0.1 - M_crit = 5 # Squared critical Mach number (above which density is modified) dim = 3 # define dimension and mesh size @@ -62,7 +61,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'te', dbc=True) + pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'te', dbc=True) tms['pre'].stop() # solve problem diff --git a/flow/tests/cylinder3.py b/flow/tests/cylinder3.py index 2a4b36c508ef55f9e2aebde61fee18028c44ea84..fe8328b07dd8440a1edb25a781af3723aa3b711b 100644 --- a/flow/tests/cylinder3.py +++ b/flow/tests/cylinder3.py @@ -41,7 +41,6 @@ def main(): alpha = 3*math.pi/180 U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1 M_inf = 0.1 - M_crit = 5 # Squared critical Mach number (above which density is modified) dim = 3 # define dimension and mesh size @@ -65,7 +64,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, dirichlet, wake, bnd,_ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'teTip', dbc=True) + pbl, dirichlet, wake, bnd,_ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'teTip', dbc=True) tms['pre'].stop() # solve problem diff --git a/flow/tests/lift.py b/flow/tests/lift.py index a5757ec99e360e07d019d7bec543e87e9b521490..52cce3bf71a097c80eab32345a439fb2c3ddec57 100644 --- a/flow/tests/lift.py +++ b/flow/tests/lift.py @@ -41,7 +41,6 @@ def main(): # define flow variables alpha = 2*math.pi/180 M_inf = 0.7 - M_crit = 5 # Squared critical Mach number (above which density is modified) c_ref = 1 dim = 2 @@ -54,7 +53,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') tms['pre'].stop() # solve problem diff --git a/flow/tests/lift3.py b/flow/tests/lift3.py index da872d0f38d0ff14e575d9d0e6b214c3c25ff771..cd05acaf71b7eee065c3b0d515efaacf582a47d8 100644 --- a/flow/tests/lift3.py +++ b/flow/tests/lift3.py @@ -40,7 +40,6 @@ def main(): # define flow variables alpha = 3*math.pi/180 M_inf = 0.3 - M_crit = 5 # Squared critical Mach number (above which density is modified) dim = 3 # define dimension and mesh size @@ -62,7 +61,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') tms['pre'].stop() # solve problem diff --git a/flow/tests/meshDef.py b/flow/tests/meshDef.py index d258e5a6d57463e480382e8ff2ac17601d16879d..f6fdba1362dd6c34064d24f921b54db79298e0a8 100644 --- a/flow/tests/meshDef.py +++ b/flow/tests/meshDef.py @@ -41,7 +41,6 @@ def main(): # define flow variables alpha0 = 0 # must be zero for this testfile M_inf = 0.3 - M_crit = 5 # Squared critical Mach number (above which density is modified) c_ref = 1 dim = 2 # parameters for mesh deformation @@ -64,9 +63,9 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, M_crit, c_ref, c_ref, xc+dx, dy, 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, c_ref, c_ref, xc+dx, dy, 0., 'airfoil', te = 'te') # set the refrence problem and add medium, initial/boundary conditions - pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, M_crit, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te') + pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te') tms['pre'].stop() # solver problem for 0°AOA diff --git a/flow/tests/meshDef3.py b/flow/tests/meshDef3.py index b1e411914db5c6bf6430bc53b853ac8a8fbab30e..11459bb6ddd8416ae8d75d12c8c3b92deb0877c8 100644 --- a/flow/tests/meshDef3.py +++ b/flow/tests/meshDef3.py @@ -40,7 +40,6 @@ def main(): # define flow variables alpha0 = 0 # must be zero for this testfile M_inf = 0.3 - M_crit = 5 # Squared critical Mach number (above which density is modified) dim = 3 # define dimension and mesh size @@ -69,7 +68,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') tms['pre'].stop() # solver problem for 0°AOA diff --git a/flow/tests/nonlift.py b/flow/tests/nonlift.py index ef6c22cc69e265943b129c1c30cb17aff91ba821..15e2bbf25532167083f6b9bdd67266397342c14b 100644 --- a/flow/tests/nonlift.py +++ b/flow/tests/nonlift.py @@ -42,7 +42,6 @@ def main(): # define flow variables alpha = 0*math.pi/180 # must be zero for this testfile M_inf = 0.8 - M_crit = 5 # Squared critical Mach number (above which density is modified) c_ref = 1 dim = 2 @@ -55,7 +54,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') tms['pre'].stop() # solve problem diff --git a/flow/validation/agard.py b/flow/validation/agard.py index edc13ba39e6505252b181d3a872c2397b8f68634..8df7bea8801347ec4a66494b2f6df31a0c9777da 100644 --- a/flow/validation/agard.py +++ b/flow/validation/agard.py @@ -70,7 +70,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') tms['pre'].stop() # solve problem @@ -100,6 +100,8 @@ def main(): tests = CTests() tests.add(CTest('CL', solver.Cl, 0.062, 1e-1)) tests.add(CTest('CD', solver.Cd, 0.0006, 1e-1)) + tests.add(CTest('CS', solver.Cs, 0.0017, 1e-1)) + tests.add(CTest('CM', solver.Cm, -0.064, 1e-1)) tests.run() # eof diff --git a/flow/validation/onera.py b/flow/validation/onera.py index ab71b2ba8b704dc8800d27ece198a7ecc2f47690..8d742f02cfee0daf2b065199c11a577ec20613e3 100644 --- a/flow/validation/onera.py +++ b/flow/validation/onera.py @@ -43,7 +43,6 @@ def main(): # define flow variables alpha = 3.06*np.pi/180 M_inf = 0.839 - M_crit = 5 # Squared critical Mach number (above which density is modified) dim = 3 # define dimension and mesh size @@ -68,7 +67,7 @@ def main(): # set the problem and add medium, initial/boundary conditions tms['pre'].start() - pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') tms['pre'].stop() # solve problem diff --git a/fwk/src/fwk.h b/fwk/src/fwk.h index dff0dee4d47dbe83349d0f258a3f1228a65c0ba7..85c99132530682b33a421d1812b568b5667f21d9 100644 --- a/fwk/src/fwk.h +++ b/fwk/src/fwk.h @@ -48,6 +48,7 @@ namespace fwk class wObject; class Extract; +class Observer; class HMSTime; class Time; diff --git a/fwk/src/wObserver.cpp b/fwk/src/wObserver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0162b4a378aaf7b7023869e4497083dccffc4339 --- /dev/null +++ b/fwk/src/wObserver.cpp @@ -0,0 +1,25 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "wObserver.h" +#include <stdexcept> + +using namespace fwk; + +void Observer::update(double data) +{ + throw std::runtime_error("fwk::Observer::update() not implemented!\n"); +} \ No newline at end of file diff --git a/fwk/src/wObserver.h b/fwk/src/wObserver.h new file mode 100644 index 0000000000000000000000000000000000000000..cabc058cc05eafab878f55a4f514490c8613e12d --- /dev/null +++ b/fwk/src/wObserver.h @@ -0,0 +1,39 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef WOBSERVER_H +#define WOBSERVER_H + +#include "fwk.h" + +namespace fwk +{ + +/** + * @brief Base class that obsverves another class, and that updates accordingly + */ + +class FWK_API Observer +{ +public: + Observer() {} + virtual ~Observer() {} + virtual void update(double data); +}; + +} // namespace fwk + +#endif //WOBSERVER_H diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp index 35bab22c0acd6493f3769b0cfbcc82630ed702c8..7bca36520e0f46825597482e9c8043974721e470 100644 --- a/heat/src/wSolver.cpp +++ b/heat/src/wSolver.cpp @@ -71,7 +71,6 @@ Solver::~Solver() } /** - * @todo detecter/traiter les jacos negatifs pour tous les éléments * @todo planter clairement si pas de Medium * @todo calculer la taille du systeme en fct du nombre de noeuds des physical groups utilisés */ diff --git a/tbox/src/wFct0.cpp b/tbox/src/wFct0.cpp index b0a7b0dc334be0a68182e6ff2885ea34927e82b3..2f4618d30a486290a8a2765062bfb4dc37ad319f 100644 --- a/tbox/src/wFct0.cpp +++ b/tbox/src/wFct0.cpp @@ -66,13 +66,6 @@ void PwLf::write(std::ostream &out) const //------------------------------------------------------------------------- -double Fct0::evalD(Element const &e, std::vector<double> const &u, size_t k) const -{ - throw std::runtime_error("Fct0::evalD not implemented\n"); -} - -//------------------------------------------------------------------------- - void Fct0C::write(std::ostream &out) const { out << "Fct0C=" << v; diff --git a/tbox/src/wFct0.h b/tbox/src/wFct0.h index 81da6cd4156c4ce04b8660aa55f1c087d675a06b..eb6d41661c5309808cb262d89b42cdea41a2406d 100644 --- a/tbox/src/wFct0.h +++ b/tbox/src/wFct0.h @@ -41,10 +41,9 @@ public: /** * @brief Scalar functions to be integrated over an element - * + * * eval() evaluates the function at kth Gauss point - * evalD() evaluates the derivative at the kth Gauss point - * @authors Romain Boman, Adrien Crovato + * @authors Romain Boman */ class TBOX_API Fct0 : public fwk::wObject { @@ -62,9 +61,6 @@ public: virtual double eval(Element const &e, std::vector<double> const &u, size_t k) const = 0; - virtual double evalD(Element const &e, - std::vector<double> const &u, - size_t k) const; #endif }; @@ -78,7 +74,6 @@ class TBOX_API Fct0C : public Fct0 public: Fct0C(double _v) : Fct0(), v(_v) {} - void update(double _v) { v = _v; } #ifndef SWIG virtual double eval(Element const &e, @@ -87,12 +82,6 @@ public: { return v; } - virtual double evalD(Element const &e, - std::vector<double> const &u, - size_t k) const override - { - return 0; - } virtual void write(std::ostream &out) const override; #endif }; diff --git a/tbox/src/wFct1.h b/tbox/src/wFct1.h index 439747f2bf3e4dc23a3e381466009067f20aa230..0a3c546dc6a54630f9a17075a1d39480a92a1ccc 100644 --- a/tbox/src/wFct1.h +++ b/tbox/src/wFct1.h @@ -27,7 +27,7 @@ namespace tbox /** * @brief Vector functions to be integrated over an element - * + * * eval() evaluates the function at kth Gauss point * @authors Romain Boman, Adrien Crovato */ @@ -47,9 +47,6 @@ public: virtual Eigen::Vector3d eval(Element const &e, std::vector<double> const &u, size_t k) const = 0; - virtual Eigen::Vector3d evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const = 0; #endif }; @@ -76,12 +73,6 @@ public: { return v; } - virtual Eigen::Vector3d evalD(tbox::Element const &e, - std::vector<double> const &u, - size_t k) const override - { - return Eigen::Vector3d::Zero(); - } virtual void write(std::ostream &out) const override; #endif };