sdpm/src/sdpmQuad4.h

@@ -37,6 +37,7 @@ public:
 #ifndef SWIG
     virtual void update() override;
     virtual void update(sdpmDouble beta) override;
+    virtual sdpmMatrixXd computeGradientOperator() const override;
     virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const override;
 #endif
 };

sdpm/src/sdpmSolver.cpp

@@ -34,7 +34,7 @@ Solver::Solver(Problem &pbl, int vrb) : _pbl(pbl), _vrb(vrb), _cl(0), _cd(0), _c
     // Say hi
     std::cout << std::endl;
     std::cout << "*******************************************************************************" << std::endl;
-    std::cout << "** Hi! My name is SDPM v1.1-2401                                             **" << std::endl;
+    std::cout << "** Hi! My name is SDPM v1.2-2412                                             **" << std::endl;
     std::cout << "** ULiege 2023-2024                                                          **" << std::endl;
     std::cout << "*******************************************************************************" << std::endl;
 
@@ -558,7 +558,7 @@ void Solver::post(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdpmVectorX
             sdpmComplex phit = -sdpmComplex(0., 1.) * omega * emu(_rows[elems[i]]);
             sdpmComplex phix = u1[i](0) - (uc0[i](0) + uc1[i](0) * sdpmComplex(0., 1.) * omega);
             sdpmComplex phix0(_u[elems[i]->getId() - 1](0) - ui(0), 0.);
-            sdpmComplex cp1 = -sdpmComplex(2., 0.) * phit + sdpmComplex(mi * mi, 0.) * (sdpmComplex(2., 0.) * phix0 * phix + phix0 * phit);
+            sdpmComplex cp1 = -sdpmComplex(2., 0.) * phit + sdpmComplex(2 * mi * mi, 0.) * (phix0 * phix + phix0 * phit);
             for (size_t j = 0; j < 3; ++j)
                 cp1 -= sdpmComplex(2 * _u[elems[i]->getId() - 1](j), 0.) * u1[i](j);
             // cp
@@ -616,7 +616,9 @@ void Solver::computeGAF(size_t ifq)
     // Get sizes
     size_t n = _pbl.getNodes().size();
     size_t nm = _pbl.getNumModes();
-    // Get z-component of modal displacements and nodal loads
+    // Get x-component of reference center
+    sdpmDouble xRef = _pbl.getRefCtr()(0);
+    // Get z-component of displacements and nodal loads
     sdpmMatrixXd mz1 = sdpmMatrixXd::Zero(nm, n);
     sdpmMatrixXd lz1Re = sdpmMatrixXd::Zero(n, nm);
     sdpmMatrixXd lz1Im = sdpmMatrixXd::Zero(n, nm);
@@ -624,12 +626,15 @@ void Solver::computeGAF(size_t ifq)
     {
         for (auto b : _pbl.getBodies())
         {
+            sdpmDouble aAmp = b->getPitchAmplitude(imd);
+            sdpmDouble hAmp = b->getPlungeAmplitude(imd);
             mz1.row(imd) = Eigen::Map<const sdpmVectorXd>(b->getModeZ(imd).data(), n).transpose();
             std::vector<Node *> nods = b->getNodes();
             std::vector<sdpmVector3d> l1Re = b->getUnsteadyLoadsReal(imd, ifq);
             std::vector<sdpmVector3d> l1Im = b->getUnsteadyLoadsImag(imd, ifq);
             for (size_t i = 0; i < nods.size(); ++i)
             {
+                mz1(imd, nods[i]->getId() - 1) += aAmp * (nods[i]->getCoords()(0) - xRef) + hAmp;
                 lz1Re(nods[i]->getId() - 1, imd) = -l1Re[i](2);
                 lz1Im(nods[i]->getId() - 1, imd) = -l1Im[i](2);
             }

sdpm/src/sdpmSolver.h

@@ -31,18 +31,22 @@ namespace sdpm
  */
 class SDPM_API Solver : public sdpmObject
 {
-    friend class Adjoint; // so that Adjoint can access the solver outputs
+    friend class Adjoint;  // so that Adjoint can access the solver outputs
+    friend class Gradient; // so that Gradient can access the solver outputs
 
+protected:
     // Problem
     Problem &_pbl; ///< problem definition
     int _vrb;      ///< verbosity level
     // AICs
     size_t _nP;                     ///< number of body panels
-    std::map<Element *, int> _rows; ///< map linking an element to a matrix cell
+    std::map<Element *, int> _rows; ///< map linking an element to a matrix entry
     sdpmMatrixXd _A0;               ///< zero-frequency body-body doublet matrix
     sdpmMatrixXd _B0;               ///< zero-frequency body-body source matrix
     std::vector<sdpmMatrixXcd> _A;  ///< unsteady body-body doublet matrix
     std::vector<sdpmMatrixXcd> _B;  ///< unsteady body-body source matrix
+
+private:
     // Solution
     std::vector<sdpmDouble> _tau;                              ///< source strength
     std::vector<sdpmDouble> _mu;                               ///< doublet strength (perturbation potential)
@@ -71,7 +75,7 @@ class SDPM_API Solver : public sdpmObject
 
 public:
     Solver(Problem &pbl, int vrb = 1);
-    ~Solver();
+    virtual ~Solver();
 
     std::vector<sdpmDouble> const &getPressure() const { return _cp; }
     std::vector<sdpmDouble> const &getUnsteadyPressureReal(size_t imd, size_t ifq) const { return _cp1Re[imd][ifq]; }
@@ -92,14 +96,14 @@ public:
     std::vector<std::vector<sdpmDouble>> getGafMatrixIm(size_t ifq) const;
 
     void update();
-    void run();
+    virtual void run();
     void save(GmshExport const &writer, std::string const &suffix = "");
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
 #endif
 
-private:
+protected:
     void post(sdpmVectorXd const &etau, sdpmVectorXd const &emu);
     void post(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu);
     void computeGAF(size_t ifq);

sdpm/src/sdpmSolverTransonic.cpp

+/*
+ * Copyright 2023 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 "sdpmSolverTransonic.h"
+#include "sdpmProblem.h"
+#include "sdpmBody.h"
+#include "sdpm_ad_interface.h"
+#include <iostream>
+
+using namespace sdpm;
+
+SolverTransonic::SolverTransonic(Problem &pbl, int vrb) : Solver(pbl, vrb), _dCorr(sdpmVectorXd::Ones(_nP))
+{
+    // Build node row id
+    int i = 0;
+    for (auto body : _pbl.getBodies())
+        for (auto n : body->getNodes())
+            if (_rowsN.insert(std::pair<Node *, int>(n, i)).second == true)
+                ++i;
+}
+
+/**
+ * @brief Solve the steady and unsteady potential equations corrected for transonic effects
+ */
+void SolverTransonic::run()
+{
+    std::cout << "Computing transonic correction...";
+    computeCorrection();
+    std::cout << "done." << std::endl;
+
+    std::cout << "Running transonic SDPM solver..." << std::flush;
+    if (_vrb > 0)
+        std::cout << "\n-- Mean flow\n"
+                  << "computing sources... " << std::flush;
+    sdpmDouble beta = _pbl.getCompressibility();
+    sdpmVector3d ui = _pbl.getDragDir();
+    sdpmVectorXd etauxy(_nP);
+    sdpmVectorXd etauz(_nP);
+    for (auto body : _pbl.getBodies())
+    {
+        for (auto e : body->getElements())
+        {
+            etauxy(_rows[e]) = ui[0] * e->getCompressibleNormal()[0] / beta + ui[1] * e->getCompressibleNormal()[1];
+            etauz(_rows[e]) = ui[2] * e->getCompressibleNormal()[2];
+        }
+    }
+    if (_vrb > 0)
+        std::cout << "done.\n"
+                  << "solving for doublets... " << std::flush;
+    sdpmVectorXd rhsxy = _B0 * etauxy;
+    sdpmVectorXd rhsz = _B0 * etauz;
+    sdpmVectorXd emuxy = solve(_A0, rhsxy);
+    sdpmVectorXd emuz = solve(_A0, rhsz);
+    if (_vrb > 0)
+        std::cout << "done.\n"
+                  << "computing flow and loads on bodies... " << std::flush;
+    post(etauxy + etauz, emuxy + emuz.cwiseProduct(_dCorr));
+    if (_vrb > 0)
+        std::cout << "done." << std::endl;
+
+    for (size_t imd = 0; imd < _pbl.getNumModes(); ++imd)
+    {
+        for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
+        {
+            if (_vrb > 0)
+                std::cout << "-- Mode #" << imd << ", frequency #" << ifq << "\n"
+                          << "computing sources... " << std::flush;
+            sdpmVectorXcd etau1(_nP);
+            sdpmVectorXcd etau1xy(_nP);
+            sdpmVectorXcd etau1z(_nP);
+            for (auto body : _pbl.getBodies())
+            {
+                std::vector<sdpmVector3d> uc0 = body->getSteadyVelocity(imd);
+                std::vector<sdpmVector3d> uc1 = body->getHarmonicVelocity(imd);
+                std::vector<Element *> elems = body->getElements();
+                for (size_t i = 0; i < elems.size(); ++i)
+                    etau1(_rows[elems[i]]) = (uc0[i] + uc1[i] * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]).conjugate().dot(elems[i]->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
+                for (size_t i = 0; i < elems.size(); ++i)
+                {
+                    etau1xy(_rows[elems[i]]) = (uc0[i](0) + uc1[i](0) * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]) * elems[i]->getCompressibleNormal()(0) / beta + (uc0[i](1) + uc1[i](1) * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]) * elems[i]->getCompressibleNormal()(1);
+                    etau1z(_rows[elems[i]]) = (uc0[i](2) + uc1[i](2) * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]) * elems[i]->getCompressibleNormal()(2);
+                }
+            }
+            if (_vrb > 0)
+                std::cout << "done.\n"
+                          << "solving for doublets... " << std::flush;
+            // Split complex set of equations into real set (twice the size)
+            sdpmMatrixXd A(2 * _nP, 2 * _nP);
+            A.block(0, 0, _nP, _nP) = _A[ifq].real();
+            A.block(_nP, _nP, _nP, _nP) = _A[ifq].real();
+            A.block(_nP, 0, _nP, _nP) = _A[ifq].imag();
+            A.block(0, _nP, _nP, _nP) = -_A[ifq].imag();
+            sdpmVectorXd rhs1xy(2 * _nP);
+            sdpmVectorXd rhs1z(2 * _nP);
+            rhs1xy.block(0, 0, _nP, 1) = (_B[ifq] * etau1xy).real();
+            rhs1xy.block(_nP, 0, _nP, 1) = (_B[ifq] * etau1xy).imag();
+            rhs1z.block(0, 0, _nP, 1) = (_B[ifq] * etau1z).real();
+            rhs1z.block(_nP, 0, _nP, 1) = (_B[ifq] * etau1z).imag();
+            sdpmVectorXd muxy = solve(A, rhs1xy);
+            sdpmVectorXd muz = solve(A, rhs1z);
+            sdpmVectorXcd emu1xy = muxy.block(0, 0, _nP, 1) + sdpmComplex(0, 1) * muxy.block(_nP, 0, _nP, 1);
+            sdpmVectorXcd emu1z = muz.block(0, 0, _nP, 1) + sdpmComplex(0, 1) * muz.block(_nP, 0, _nP, 1);
+            if (_vrb > 0)
+                std::cout << "done.\n"
+                          << "computing flow and loads on bodies... " << std::flush;
+            sdpmVectorXcd mu1c = emu1xy + _dCorr.cwiseProduct(emu1z);
+            post(imd, ifq, etau1, mu1c);
+            if (_vrb > 0)
+                std::cout << "done." << std::endl;
+        }
+    }
+
+    // GAF
+    if (_vrb > 0)
+        std::cout << "-- Generalized Aerodynamic Force matrices" << std::endl;
+    for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
+    {
+        if (_vrb > 0)
+            std::cout << "computing GAF at frequency #" << ifq << "... " << std::flush;
+        computeGAF(ifq);
+        if (_vrb > 0)
+            std::cout << "done." << std::endl;
+    }
+    if (_vrb == 0)
+        std::cout << " done.";
+    std::cout << std::endl;
+}
+
+/**
+ * @brief Compute the transonic correction terms
+ */
+void SolverTransonic::computeCorrection()
+{
+    // Assemble x-gradient matrix
+    // 1) Gradient matrix on element
+    sdpmMatrixXd dNx = sdpmMatrixXd::Zero(_rows.size(), _rowsN.size());
+    for (auto p : _rows)
+    {
+        Element *e = p.first;
+        std::vector<Node *> nodes = e->getNodes();
+        sdpmVectorXd dNe = e->computeGradientOperator().row(0).transpose();
+        for (size_t in = 0; in < nodes.size(); ++in)
+            dNx(_rows.at(e), _rowsN.at(nodes[in])) = dNe(in);
+    }
+    // 2) Average matrix from element to node
+    sdpmMatrixXd L = sdpmMatrixXd::Zero(_rowsN.size(), _rows.size());
+    for (auto body : _pbl.getBodies())
+    {
+        std::map<Node *, std::vector<Element *>> neMap = body->getMap();
+        for (auto n : body->getNodes())
+            for (auto e : neMap.at(n))
+                L(_rowsN.at(n), _rows.at(e)) = 1. / neMap.at(n).size();
+    }
+    sdpmMatrixXd dNxL = dNx * L;
+    // 3) Boundary condition on upper trailing edge
+    for (auto body : _pbl.getBodies())
+    {
+        std::vector<Element *> teElems;
+        std::tie(teElems, std::ignore) = body->getTrailingEdgeElements();
+        for (auto e : teElems)
+        {
+            dNxL.row(_rows.at(e)) = sdpmVectorXd::Zero(_nP).transpose();
+            dNxL(_rows.at(e), _rows.at(e)) = 1.0;
+        }
+    }
+
+    // Assemble RHS
+    // 1) Average derivative of pressure wrt angle of attack at element
+    sdpmVectorXd dCpElems = sdpmVectorXd::Zero(_nP);
+    for (auto body : _pbl.getBodies())
+    {
+        std::vector<sdpmDouble> dCpNodes = body->getTransonicPressureGrad();
+        for (auto e : body->getElements())
+            for (auto n : e->getNodes())
+                dCpElems(_rows.at(e)) += dCpNodes[n->getId() - 1] / e->getNodes().size();
+    }
+    // 2) Compute derivative of doublet wrt angle of attack
+    sdpmVectorXd nz(_nP);
+    for (auto body : _pbl.getBodies())
+        for (auto e : body->getElements())
+            nz(_rows[e]) = e->getCompressibleNormal()[2] * cos(_pbl.getAos());
+    sdpmVectorXd rhs = _B0 * nz;
+    sdpmVectorXd dMuA = solve(_A0, rhs);
+    // 3) Compute RHS with boundary condition
+    sdpmVectorXd dMuAxCorr = sdpmVectorXd::Zero(_nP);
+    for (auto body : _pbl.getBodies())
+    {
+        for (auto e : body->getElements())
+            dMuAxCorr(_rows.at(e)) = 0.5 * _pbl.getCompressibility() * dCpElems(_rows.at(e)) - nz(_rows[e]) * e->getCompressibleNormal()[0];
+        std::vector<Element *> teElems;
+        std::tie(teElems, std::ignore) = body->getTrailingEdgeElements();
+        for (auto e : teElems)
+            dMuAxCorr(_rows.at(e)) = dMuA(_rows.at(e));
+    }
+
+    // Solve for correction
+    sdpmVectorXd dMuACorr = solve(dNxL, dMuAxCorr);
+    _dCorr = dMuACorr.cwiseQuotient(dMuA);
+}
+
+void SolverTransonic::write(std::ostream &out) const
+{
+    out << "sdpm::SolverTransonic with\n"
+        << _pbl;
+}

sdpm/src/sdpmSolverTransonic.h

+/*
+ * Copyright 2023 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 SDPMSOLVERTRANSONIC_H
+#define SDPMSOLVERTRANSONIC_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+#include "sdpmSolver.h"
+#include <map>
+
+namespace sdpm
+{
+
+/**
+ * @brief Source-doublet panel solver corrected for transonic effects
+ * @authors Adrien Crovato
+ */
+class SDPM_API SolverTransonic : public Solver
+{
+    std::map<Node *, int> _rowsN; ///< map linking a node to a matrix entry
+    sdpmVectorXd _dCorr;          ///< transonic correction terms
+
+public:
+    SolverTransonic(Problem &pbl, int vrb = 1);
+
+    void run() override;
+
+#ifndef SWIG
+    virtual void write(std::ostream &out) const override;
+#endif
+
+private:
+    void computeCorrection();
+};
+
+} // namespace sdpm
+#endif // SDPMSOLVERTRANSONIC_H

sdpm/src/sdpmWake.cpp

@@ -16,6 +16,7 @@
 
 #include "sdpmWake.h"
 #include "sdpmEdge.h"
+#include "sdpmMesh.h"
 #include "sdpmNode.h"
 #include "sdpmElement.h"
 #include <vector>
@@ -35,6 +36,19 @@ Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, size_t nTe,
         teEdges.insert(ed);
     }
 
+    // Save nodes
+    std::vector<Node *> mshNodes = msh.getNodes();
+    size_t endNodW = mshNodes.size();
+    size_t nNodWake = ((elems.size() / nTe) + 1) * (nTe + 1);
+    _nodes.resize(nNodWake);
+    // first nTe+1 nodes are trailing edge nodes
+    for (size_t i = 0; i <= nTe; ++i)
+        _nodes[i] = elems[i]->getNodes()[0];
+    _nodes[nTe] = elems[nTe - 1]->getNodes()[3];
+    // other nodes are last added nodes in the mesh (other than duplicated trailing edge nodes)
+    for (size_t i = 1; i <= nNodWake - nTe - 1; ++i)
+        _nodes[nNodWake - i] = mshNodes[endNodW - i];
+
     // Find wing elements having an edge on the TE...
     // ... the elements on the pressure side still have their original (suction side) trailing edge nodes
     if (iTeMap.empty())

sdpm/src/sdpmWake.h

@@ -31,6 +31,7 @@ namespace sdpm
  */
 class SDPM_API Wake : public Group
 {
+    std::vector<Node *> _nodes;                                  ///< nodes of the wake
     std::map<Element *, std::pair<Element *, Element *>> _wwMap; ///< wing-wake map
     std::map<Element *, sdpmDouble> _lagMap;                     ///< time lag map
 
@@ -39,10 +40,11 @@ public:
     ~Wake() {}
 
 #ifndef SWIG
-    std::map<Element *, std::pair<Element *, Element *>> const &getElMap() const
+    std::vector<Node *> const &getNodes() const
     {
-        return _wwMap;
+        return _nodes;
     }
+    std::map<Element *, std::pair<Element *, Element *>> const &getElMap() const { return _wwMap; }
     std::map<Element *, sdpmDouble> const &getLagMap() const { return _lagMap; }
     virtual void write(std::ostream &out) const override;
 #endif

sdpm/tests/agard.py

@@ -76,11 +76,11 @@ def main():
     tests.add(CTest('Abs(CL1), mode 1', abs(cl1[0]), 0.076, 5e-2))
     tests.add(CTest('Ang(CL1), mode 1', cmath.phase(cl1[0]) * 180 / math.pi, 7.6, 5e-2))
     tests.add(CTest('Abs(CL1), mode 2', abs(cl1[1]), 0.114, 5e-2))
-    tests.add(CTest('Ang(CL1), mode 2', cmath.phase(cl1[1]) * 180 / math.pi, -4.7, 5e-2))
+    tests.add(CTest('Ang(CL1), mode 2', cmath.phase(cl1[1]) * 180 / math.pi, -4.8, 5e-2))
     tests.add(CTest('Abs(CL1), mode 3', abs(cl1[2]), 0.099, 5e-2))
     tests.add(CTest('Ang(CL1), mode 3', cmath.phase(cl1[2]) * 180 / math.pi, 177.0, 5e-2))
     tests.add(CTest('Abs(CL1), mode 4', abs(cl1[3]), 0.015, 5e-2))
-    tests.add(CTest('Ang(CL1), mode 4', cmath.phase(cl1[3]) * 180 / math.pi, 176.4, 5e-2))
+    tests.add(CTest('Ang(CL1), mode 4', cmath.phase(cl1[3]) * 180 / math.pi, 176.1, 5e-2))
     tests.run()
 
     # eof

sdpm/tests/agard_grad.py

+# -*- coding: utf-8 -*-
+
+# Copyright 2023 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.
+
+# AGARD 445.6 wing (gradient)
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath, interpolate_modes
+from sdpm.testing import *
+import numpy as np
+
+def main():
+    # geometry parameters
+    wnc = 20 # number of chordwise panels
+    pars = {'nC': wnc, 'nS': 10, 'pC': 1.1, 'pS': 1.0}
+    sref = 0.353 # reference area
+    crd = 0.559 # root chord
+    # flow and time parameters
+    minf = 0.9 # freestream Mach number
+    aoa = 1.0 * np.pi / 180 # angle of attack
+    kred = 0.1 # reduced frequency
+    n_modes = 4 # number of modes
+
+    # start timer
+    tms = sdpm.Timers()
+    tms['total'].start()
+    # mesh
+    tms['msh'].start()
+    msh = MeshLoader(build_fpath('models/agard445.geo', __file__, 1)).run(pars)
+    tms['msh'].stop()
+    # problem
+    tms['pbl'].start()
+    pbl = sdpm.Problem(msh)
+    pbl.setGeometry(sref, crd, 0., 0., 0., True)
+    pbl.setFreestream(aoa, 0., minf)
+    pbl.setUnsteady(n_modes, [kred])
+    wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, n_modes, 1)
+    x_modes, amp_modes = interpolate_modes(build_fpath('models/agard445_modes.csv', __file__, 1), wing.getNodes())
+    for i in range(n_modes):
+        wing.addMotion()
+        wing.setFlexibleMotion(i, 1 / amp_modes[i], x_modes[:, 3*i:3*i+3])
+    pbl.addBody(wing)
+    tms['pbl'].stop()
+    # solver
+    tms['sol'].start()
+    sol = sdpm.Solver(pbl)
+    sol.update()
+    sol.run()
+    sol.save(sdpm.GmshExport(msh))
+    tms['sol'].stop()
+    # gradient evaluation
+    tms['grd'].start()
+    grd = sdpm.Gradient(sol)
+    grd.run()
+    tms['grd'].stop()
+    # end timer
+    tms['total'].stop()
+    print(tms)
+
+    # recover GAF matrix from partial gradients (NB: since GAF is quadratic in DOF, result must be divided by 2)
+    q_re = np.array(sol.getGafMatrixRe(0), dtype=float)
+    q_im = np.array(sol.getGafMatrixIm(0), dtype=float)
+    dq_re = np.zeros((n_modes, n_modes), dtype=float)
+    dq_im = np.zeros((n_modes, n_modes), dtype=float)
+    for i in range(n_modes):
+        for j in range(n_modes):
+            dq_re[i,j] = np.array(grd.computePartialsGafRe(0, j, i, j, 'dz'), dtype=float).dot(x_modes[:, 3*j+0] / amp_modes[j]) \
+                       + np.array(grd.computePartialsGafRe(0, j, i, j, 'rx'), dtype=float).dot(x_modes[:, 3*j+1] / amp_modes[j]) \
+                       + np.array(grd.computePartialsGafRe(0, j, i, j, 'ry'), dtype=float).dot(x_modes[:, 3*j+2] / amp_modes[j])
+            dq_im[i,j] = np.array(grd.computePartialsGafIm(0, j, i, j, 'dz'), dtype=float).dot(x_modes[:, 3*j+0] / amp_modes[j]) \
+                       + np.array(grd.computePartialsGafIm(0, j, i, j, 'rx'), dtype=float).dot(x_modes[:, 3*j+1] / amp_modes[j]) \
+                       + np.array(grd.computePartialsGafIm(0, j, i, j, 'ry'), dtype=float).dot(x_modes[:, 3*j+2] / amp_modes[j])
+            if i != j:
+                dq_re[i,j] += np.array(grd.computePartialsGafRe(0, i, i, j, 'dz'), dtype=float).dot(x_modes[:, 3*i] / amp_modes[i])
+                dq_im[i,j] += np.array(grd.computePartialsGafIm(0, i, i, j, 'dz'), dtype=float).dot(x_modes[:, 3*i] / amp_modes[i])
+    dq_re *= 0.5
+    dq_im *= 0.5
+
+    # test
+    tests = CTests()
+    tests.add(CTest('Re(Q)',np.linalg.norm(q_re - dq_re), 0.0, 1e-10, forceabs=True))
+    tests.add(CTest('Im(Q)',np.linalg.norm(q_im - dq_im), 0.0, 1e-10, forceabs=True))
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()

sdpm/tests/lann.py

@@ -52,7 +52,7 @@ def main():
     pbl.setUnsteady(1, [kred])
     wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, 1, 1)
     wing.addMotion()
-    wing.setRigidMotion(0, amp, 0., xf/math.sqrt(1-minf*minf), 0.)
+    wing.setRigidMotion(0, amp, 0.)
     pbl.addBody(wing)
     tms['pbl'].stop()
     # solver
@@ -76,14 +76,14 @@ def main():
     tests.add(CTest('CD0', sol.getDragCoeff(),  0.0020, 5e-2)) # MATLAB: 0.0031
     tests.add(CTest('CS0', sol.getSideforceCoeff(), 0.006, 5e-2))
     tests.add(CTest('CM0', sol.getMomentCoeff(), -0.079, 5e-2))
-    tests.add(CTest('Abs(CL1)', abs(cl1) / amp / math.pi, 1.61, 5e-2)) # MATLAB: 1.69
-    tests.add(CTest('Ang(CL1)', cmath.phase(cl1) * 180 / math.pi, 0.45, 5e-2)) # MATLAB: 0.64
-    tests.add(CTest('Abs(CD1)', abs(cd1) / amp / math.pi, 0.036, 5e-2)) # MATLAB: 0.037
-    tests.add(CTest('Ang(CD1)', cmath.phase(cd1) * 180 / math.pi, 16.4, 5e-2)) # MATLAB: 19.9
+    tests.add(CTest('Abs(CL1)', abs(cl1) / amp / math.pi, 1.60, 5e-2))
+    tests.add(CTest('Ang(CL1)', cmath.phase(cl1) * 180 / math.pi, -1.8, 5e-2))
+    tests.add(CTest('Abs(CD1)', abs(cd1) / amp / math.pi, 0.036, 5e-2))
+    tests.add(CTest('Ang(CD1)', cmath.phase(cd1) * 180 / math.pi, 16.5, 5e-2))
     tests.add(CTest('Abs(CS1)', abs(cs1) / amp / math.pi, 0.040, 5e-2))
-    tests.add(CTest('Ang(CS1)', cmath.phase(cs1) * 180 / math.pi,  -9.7, 5e-2))
-    tests.add(CTest('Abs(CM1)', abs(cm1) / amp / math.pi, 0.34, 5e-2))
-    tests.add(CTest('Ang(CM1)', cmath.phase(cm1) * 180 / math.pi, -170.4, 5e-2))
+    tests.add(CTest('Ang(CS1)', cmath.phase(cs1) * 180 / math.pi,  -10.3, 5e-2))
+    tests.add(CTest('Abs(CM1)', abs(cm1) / amp / math.pi, 0.33, 5e-2))
+    tests.add(CTest('Ang(CM1)', cmath.phase(cm1) * 180 / math.pi, -174.7, 5e-2))
     tests.run()
 
     # eof

sdpm/tests/lann_ad.py

@@ -52,7 +52,7 @@ def main():
     pbl.setUnsteady(1, [kred])
     wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, 1, 1)
     wing.addMotion()
-    wing.setRigidMotion(0, amp, 0., xf/math.sqrt(1-minf*minf), 0.)
+    wing.setRigidMotion(0, amp, 0.)
     pbl.addBody(wing)
     tms['pbl'].stop()
     # solver
@@ -75,14 +75,14 @@ def main():
 
     # test
     tests = CTests()
-    tests.add(CTest('dRe(CL1)/dAlpha', d_cl1_re['aoa'][0], 1.4e-3, 2e-4, forceabs=True))
-    tests.add(CTest('dRe(CL1)/dOmega', d_cl1_re['omega'][0], -1.3e-3, 2e-4, forceabs=True))
+    tests.add(CTest('dRe(CL1)/dAlpha', d_cl1_re['aoa'][0], 1.3e-3, 2e-4, forceabs=True))
+    tests.add(CTest('dRe(CL1)/dOmega', d_cl1_re['omega'][0], -1.4e-3, 2e-4, forceabs=True))
     tests.add(CTest('dIm(CL1)/dAlpha', d_cl1_im['aoa'][0], -1.4e-4, 2e-5, forceabs=True))
-    tests.add(CTest('dIm(CL1)/dOmega', d_cl1_im['omega'][0], 1.1e-3, 2e-4, forceabs=True))
-    tests.add(CTest('dRe(CD1)/dAlpha', d_cd1_re['aoa'][0], 7.9e-3, 2e-4, forceabs=True))
-    tests.add(CTest('dRe(CD1)/dOmega', d_cd1_re['omega'][0], 3.3e-5, 2e-6, forceabs=True))
-    tests.add(CTest('dIm(CD1)/dAlpha', d_cd1_im['aoa'][0], 1.7e-3, 2e-4, forceabs=True))
-    tests.add(CTest('dIm(CD1)/dOmega', d_cd1_im['omega'][0], 6.2e-5, 2e-6, forceabs=True))
+    tests.add(CTest('dIm(CL1)/dOmega', d_cl1_im['omega'][0], 6.4e-4, 2e-4, forceabs=True))
+    tests.add(CTest('dRe(CD1)/dAlpha', d_cd1_re['aoa'][0], 8.0e-3, 2e-4, forceabs=True))
+    tests.add(CTest('dRe(CD1)/dOmega', d_cd1_re['omega'][0], 3.7e-5, 2e-6, forceabs=True))
+    tests.add(CTest('dIm(CD1)/dAlpha', d_cd1_im['aoa'][0], 1.6e-3, 2e-4, forceabs=True))
+    tests.add(CTest('dIm(CD1)/dOmega', d_cd1_im['omega'][0], 6.4e-5, 2e-6, forceabs=True))
     tests.run()
 
     # eof

sdpm/tests/lann_grad.py

+# -*- coding: utf-8 -*-
+
+# Copyright 2023 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.
+
+# LANN wing (gradient)
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath
+from sdpm.testing import *
+import numpy as np
+
+def main():
+    # geometry parameters
+    wnc = 20 # number of chordwise panels
+    pars = {'nC': wnc, 'bC': 0.25, 'rS': 1}
+    sref = 0.253 # reference area
+    crd = 0.361 # root chord
+    xf = 0.224 # x-coordinate of rotation center
+    # flow and time parameters
+    minf = 0.621 # freestream Mach number
+    aoa = 0.59 * np.pi / 180
+    kred = 0.133 # reduced frequency
+    a_amp = 0.5 * 0.25 * np.pi / 180 # semi pitching amplitude
+    h_amp = 0.1 # plunging amplitude
+
+    # start timer
+    tms = sdpm.Timers()
+    tms['total'].start()
+    # mesh
+    tms['msh'].start()
+    msh = MeshLoader(build_fpath('models/lann.geo', __file__, 1)).run(pars)
+    tms['msh'].stop()
+    # problem
+    tms['pbl'].start()
+    pbl = sdpm.Problem(msh)
+    pbl.setGeometry(sref, crd, xf, 0., 0., True)
+    pbl.setFreestream(aoa, 0., minf)
+    pbl.setUnsteady(1, [kred])
+    wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, 1, 1)
+    wing.addMotion()
+    wing.setRigidMotion(0, a_amp, h_amp)
+    pbl.addBody(wing)
+    tms['pbl'].stop()
+    # solver
+    tms['sol'].start()
+    sol = sdpm.Solver(pbl)
+    sol.update()
+    sol.run()
+    sol.save(sdpm.GmshExport(msh))
+    tms['sol'].stop()
+    # gradient evaluation
+    tms['grd'].start()
+    grd = sdpm.Gradient(sol)
+    grd.run()
+    tms['grd'].stop()
+    # end timer
+    tms['total'].stop()
+    print(tms)
+
+    # recover GAF matrix from partial gradients (NB: since GAF is quadratic in DOF, result must be divided by 2)
+    q_re = float(sol.getGafMatrixRe(0)[0][0])
+    q_im = float(sol.getGafMatrixIm(0)[0][0])
+    dq_re = 0.5 * (float(grd.computePartialsGafRe(0, 0, 0, 0, 'a')[0]) * a_amp + float(grd.computePartialsGafRe(0, 0, 0, 0, 'h')[0]) * h_amp)
+    dq_im = 0.5 * (float(grd.computePartialsGafIm(0, 0, 0, 0, 'a')[0]) * a_amp + float(grd.computePartialsGafIm(0, 0, 0, 0, 'h')[0]) * h_amp)
+
+    # test
+    tests = CTests()
+    tests.add(CTest('Re(Q)', q_re - dq_re, 0.0, 1e-10, forceabs=True))
+    tests.add(CTest('Im(Q)', q_im - dq_im, 0.0, 1e-10, forceabs=True))
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()

sdpm/tests/lann_xsonic.py

+# -*- coding: utf-8 -*-
+
+# Copyright 2023 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.
+
+# LANN wing (transonic)
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath, interpolate_pressure
+from sdpm.testing import *
+import math
+import cmath
+
+def main():
+    # geometry parameters
+    wnc = 20 # number of chordwise panels
+    pars = {'nC': wnc, 'bC': 0.25, 'rS': 1}
+    sref = 0.253 # reference area
+    crd = 0.361 # root chord
+    xf = 0.224 # x-coordinate of rotation center
+    # flow and time parameters
+    minf = 0.771 # freestream Mach number
+    aoa = 2.6 * math.pi / 180 # angle of attack
+    kred = 0.108 # reduced frequency
+    amp = 0.5 * 0.25 * math.pi / 180 # semi pitching amplitude
+
+    # start timer
+    tms = sdpm.Timers()
+    tms['total'].start()
+    # mesh
+    tms['msh'].start()
+    msh = MeshLoader(build_fpath('models/lann.geo', __file__, 1)).run(pars)
+    tms['msh'].stop()
+    # problem
+    tms['pbl'].start()
+    pbl = sdpm.Problem(msh)
+    pbl.setGeometry(sref, crd, xf, 0., 0., True)
+    pbl.setFreestream(aoa, 0., minf)
+    pbl.setUnsteady(1, [kred])
+    wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, 1, 1)
+    wing.addMotion()
+    wing.setRigidMotion(0, amp, 0.)
+    dcp = interpolate_pressure(build_fpath('models/lann_dcp_rans.csv', __file__, 1), wing.getNodes())
+    wing.setTransonicPressureGrad(dcp[:,0])
+    pbl.addBody(wing)
+    tms['pbl'].stop()
+    # solver
+    tms['sol'].start()
+    sol = sdpm.SolverTransonic(pbl)
+    sol.update()
+    sol.run()
+    sol.save(sdpm.GmshExport(msh))
+    tms['sol'].stop()
+    # end timer
+    tms['total'].stop()
+    print(tms)
+
+    # test
+    cl1 = sol.getUnsteadyLiftCoeffReal(0, 0) + 1j * sol.getUnsteadyLiftCoeffImag(0, 0)
+    cd1 = sol.getUnsteadyDragCoeffReal(0, 0) + 1j * sol.getUnsteadyDragCoeffImag(0, 0)
+    cs1 = sol.getUnsteadySideforceCoeffReal(0, 0) + 1j * sol.getUnsteadySideforceCoeffImag(0, 0)
+    cm1 = sol.getUnsteadyMomentCoeffReal(0, 0) + 1j * sol.getUnsteadyMomentCoeffImag(0, 0)
+    tests = CTests()
+    tests.add(CTest('CL0', sol.getLiftCoeff(), 0.58, 5e-2))
+    tests.add(CTest('CD0', sol.getDragCoeff(),  0.0127, 5e-2))
+    tests.add(CTest('CS0', sol.getSideforceCoeff(), 0.015, 5e-2))
+    tests.add(CTest('CM0', sol.getMomentCoeff(), -0.148, 5e-2))
+    tests.add(CTest('Abs(CL1)', abs(cl1) / amp / math.pi, 2.42, 5e-2))
+    tests.add(CTest('Ang(CL1)', cmath.phase(cl1) * 180 / math.pi, -6.5, 5e-2))
+    tests.add(CTest('Abs(CD1)', abs(cd1) / amp / math.pi, 0.192, 5e-2))
+    tests.add(CTest('Ang(CD1)', cmath.phase(cd1) * 180 / math.pi, 2.8, 5e-2))
+    tests.add(CTest('Abs(CS1)', abs(cs1) / amp / math.pi, 0.072, 5e-2))
+    tests.add(CTest('Ang(CS1)', cmath.phase(cs1) * 180 / math.pi, -13.5, 5e-2))
+    tests.add(CTest('Abs(CM1)', abs(cm1) / amp / math.pi, 0.50, 5e-2))
+    tests.add(CTest('Ang(CM1)', cmath.phase(cm1) * 180 / math.pi, -179.1, 5e-2))
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()

sdpm/utils.py

@@ -91,9 +91,12 @@ def create_slices(body, sol, fname, ys):
     for i in range(len(ys)):
         # get data on cut section
         coords, connct, cvalues = extract_data(cut_grid(grid, [0., ys[i], 0.]), snames)
-        # normalize x-coordinate by local chord
-        xc = np.zeros((coords.shape[0], 1))
-        xc[:,0] = (coords[:,0] - min(coords[:,0])) / (max(coords[:,0]) - min(coords[:,0]))
+        # normalize x and z-coordinates by local chord
+        xc = np.zeros((coords.shape[0], 2))
+        ile = np.argmin(coords[:,0]) # LE index
+        crd = max(coords[:,0]) - min(coords[:,0]) # chord
+        xc[:,0] = (coords[:,0] - coords[ile,0]) / crd
+        xc[:,1] = (coords[:,2] - coords[ile,2]) / crd
         # average value at points from values at cell center
         nvalues = np.zeros((coords.shape[0], len(sol.keys())))
         for j in range(connct.shape[0]):
@@ -101,16 +104,16 @@ def create_slices(body, sol, fname, ys):
             for k in range(len(pids)):
                 for l in range(len(sol.keys())):
                     nvalues[connct[j][k], l] += cvalues[snames[l]][j] / len(pids)
-        # create dataset holding x-coordinates and values and sort it
+        # create dataset holding coordinates and values and sort it
         xc_vals = np.hstack((xc, nvalues))
         sort(connct, xc_vals)
         xc_vals = np.vstack((xc_vals, xc_vals[0,:]))
         # write to file
-        hdr = 'x/c'
+        hdr = f'y = {ys[i]}, c = {crd}\n'
+        hdr += '{:>9s}, {:>10s}'.format('x/c', 'z/c')
         for name in snames:
-            hdr += f', {name}'
-        hdr += f'\ny = {ys[i]}'
-        np.savetxt(fname+f'_slice_{i}.dat', xc_vals, fmt='%1.5e', delimiter=',', header=hdr)
+            hdr += ', {:>10s}'.format(name)
+        np.savetxt(fname+f'_slice_{i}.dat', xc_vals, fmt='%+1.4e', delimiter=',', header=hdr)
 
 def interpolate_modes(fname, surf_nodes, flat=True):
     """Read modal displacements and interpolate them on surface grid
@@ -151,3 +154,32 @@ def interpolate_modes(fname, surf_nodes, flat=True):
     for i in range(n_modes):
         iamp[i] = np.linalg.norm(np.reshape(imodes[:, 3*i:3*i+3], (3 * len(surf_nodes), 1)))
     return imodes, iamp
+
+def interpolate_pressure(fname, surf_nodes):
+    """Read pressure derivative and interpolate it on surface grid
+
+    Parameters
+    ----------
+    fname : string
+        name of the CSV file containing the pressure derivative
+    surf_nodes : array
+        nodes of surface grid
+
+    Return
+    ------
+    icp : (n, 1) ndarray
+        interpolated pressure derivative (n is the number of surf_nodes)
+    """
+    import numpy as np
+    import scipy.interpolate as spip
+    # Read pressure and remove duplicates
+    cp = np.loadtxt(fname, delimiter=',', skiprows=1)
+    _, idx = np.unique(cp[:,:3], axis=0, return_index=True)
+    cp = cp[idx,:]
+    # Interpolate on surface mesh
+    xyz_surf = np.zeros((len(surf_nodes), 3))
+    for i, n in enumerate(surf_nodes):
+        for j in range(3):
+            xyz_surf[i,j] = n.getCoords()[j]
+    icp = spip.RBFInterpolator(cp[:,:3], cp[:,3:], smoothing=1e-6)(xyz_surf)
+    return icp