CMake/sdpmMacros.cmake

 # 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.
@@ -18,14 +18,17 @@ MACRO(MACRO_AddTest srcDir)
     file(GLOB tfiles RELATIVE ${srcDir} ${srcDir}/*)
     foreach(tfile ${tfiles})
         set(spath ${srcDir}/${tfile})
-        if(NOT IS_DIRECTORY ${spath} AND 
-           ${spath} MATCHES ".+\\.py$" AND 
-           NOT ${tfile} STREQUAL "__init__.py")
-            string(REPLACE "${PROJECT_SOURCE_DIR}/" "" strip ${spath}) 
-            message(STATUS "Adding test ${strip}")
-            add_test(NAME ${strip} 
-                    WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} 
-                    COMMAND ${Python3_EXECUTABLE} "${PROJECT_SOURCE_DIR}/run.py" --nogui --clean ${strip})            
+        if(NOT IS_DIRECTORY ${spath} AND ${spath} MATCHES ".+\\.py$" AND NOT ${tfile} STREQUAL "__init__.py")
+            if(${USE_CODI} MATCHES "OFF" AND ${tfile} MATCHES "._ad+\\.py$")
+                string(REPLACE "${PROJECT_SOURCE_DIR}/" "" strip ${spath})
+                message(STATUS "Skipping test ${strip}")
+            else()
+                string(REPLACE "${PROJECT_SOURCE_DIR}/" "" strip ${spath})
+                message(STATUS "Adding test ${strip}")
+                add_test(NAME ${strip}
+                        WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+                        COMMAND ${Python3_EXECUTABLE} "${PROJECT_SOURCE_DIR}/run.py" --nogui --clean ${strip})
+            endif()
         else()
             MACRO_AddTest(${srcDir}/${tfile})
         endif()

CMakeLists.txt

@@ -56,7 +56,7 @@ ELSEIF(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
 ENDIF()
 
 # -- OPTIONS
-OPTION(USE_CODI         "Compile with CoDiPack"         ON)
+OPTION(USE_CODI         "Compile with CoDiPack"         OFF)
 
 # -- DEPENDENCIES
 # Python

sdpm/_src/sdpmw.i

@@ -43,6 +43,7 @@ threads="1"
 #include "sdpmSolver.h"
 #include "sdpmSolverTransonic.h"
 #include "sdpmAdjoint.h"
+#include "sdpmGradient.h"
 
 #include "sdpmCppBuf2Py.h"
 
@@ -140,6 +141,7 @@ typedef double sdpmDouble; // needed so that SWIG does not convert sdpmDouble to
 %include "sdpmSolver.h"
 %include "sdpmSolverTransonic.h"
 %include "sdpmAdjoint.h"
+%include "sdpmGradient.h"
 
 %warnfilter(401); // Nothing known about base class 'std::basic_streambuf< char >'
 %include "sdpmCppBuf2Py.h"

sdpm/src/sdpm.h

@@ -81,6 +81,7 @@ class Builder;
 class Solver;
 class SolverTransonic;
 class Adjoint;
+class Gradient;
 } // namespace sdpm
 
 #endif // SDPM_H

sdpm/src/sdpmGradient.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 "sdpmGradient.h"
+#include "sdpmSolver.h"
+#include "sdpmProblem.h"
+#include "sdpmBody.h"
+#include <iostream>
+
+using namespace sdpm;
+
+Gradient::Gradient(Solver &sol) : _sol(sol),
+                                  _nF(_sol._pbl._omega.size()),
+                                  _nM(_sol._pbl.getNumModes())
+{
+    // Create node row identification map
+    int inod = 0;
+    for (auto body : _sol._pbl.getBodies())
+        for (auto n : body->getNodes())
+            _rowsN.insert(std::pair<Node *, int>(n, inod++));
+    _nN = _rowsN.size();
+
+    // Resize displacement and load containers
+    _dz.resize(_nM, sdpmVectorXd::Zero(_nN));
+    _clz1Re.resize(_nM, std::vector<sdpmVectorXd>(_nF, sdpmVectorXd::Zero(_nN)));
+    _clz1Im.resize(_nM, std::vector<sdpmVectorXd>(_nF, sdpmVectorXd::Zero(_nN)));
+
+    // Resize pressure gradient containers
+    _dCp1Re_a.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_a.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_rx.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_rx.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_ry.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_ry.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_aDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_aDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_hDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_hDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_rxDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_rxDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_ryDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_ryDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_dzDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_dzDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_aDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_aDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_hDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_hDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_dzDotDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_dzDotDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+
+    // Resize load gradient containers
+    _dClz1Re_rx.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Im_rx.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Re_ry.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Im_ry.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Re_dz.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Im_dz.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+}
+
+Gradient::~Gradient()
+{
+    std::cout << "~sdpm::Gradient()\n";
+}
+
+/**
+ * @brief Compute partial derivatives of pressure and loads with respect to degrees of freedom
+ */
+void Gradient::run()
+{
+    // Store loads and displacements
+    for (auto body : _sol._pbl.getBodies())
+    {
+        for (size_t imd = 0; imd < _nM; ++imd)
+        {
+            std::vector<sdpmDouble> dz = body->getModeZ(imd);
+            for (auto n : body->getNodes())
+                _dz[imd](_rowsN.at(n)) = dz[n->getId() - 1];
+            for (size_t ifq = 0; ifq < _nF; ++ifq)
+            {
+                std::vector<sdpmVector3d> clzre = body->getUnsteadyLoadsReal(imd, ifq);
+                std::vector<sdpmVector3d> clzim = body->getUnsteadyLoadsImag(imd, ifq);
+                std::vector<Node *> nods = body->getNodes();
+                for (size_t i = 0; i < nods.size(); ++i)
+                {
+                    _clz1Re[imd][ifq](_rowsN.at(nods[i])) = clzre[i](2);
+                    _clz1Im[imd][ifq](_rowsN.at(nods[i])) = clzim[i](2);
+                }
+            }
+        }
+    }
+
+    // Surface gradient operator and average matrices
+    std::vector<sdpmMatrixXd> dN;
+    dN.resize(3, sdpmMatrixXd::Zero(_sol._rows.size(), _rowsN.size()));
+    for (auto p : _sol._rows)
+    {
+        Element *e = p.first;
+        std::vector<Node *> nodes = e->getNodes();
+        sdpmMatrixXd dNe = e->computeGradientOperator();
+        for (size_t in = 0; in < nodes.size(); ++in)
+            for (size_t idim = 0; idim < 3; ++idim)
+                dN[idim](_sol._rows.at(e), _rowsN.at(nodes[in])) = dNe(idim, in);
+    }
+    sdpmMatrixXd L = sdpmMatrixXd::Zero(_rowsN.size(), _sol._rows.size()); // elements-to-node average
+    for (auto body : _sol._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), _sol._rows.at(e)) = 1. / neMap.at(n).size();
+    }
+    sdpmMatrixXd L1 = sdpmMatrixXd::Zero(_sol._rows.size(), _rowsN.size()); // nodes-to-element average
+    for (auto body : _sol._pbl.getBodies())
+        for (auto e : body->getElements())
+            for (auto n : e->getNodes())
+                L1(_sol._rows.at(e), _rowsN.at(n)) = 1. / e->getNodes().size();
+
+    // Freestream values
+    sdpmDouble beta = _sol._pbl.getCompressibility();           // compressiblity factor
+    sdpmDouble mach = _sol._pbl.getMach();                      // Mach number
+    sdpmVector3d uinf = _sol._pbl.getDragDir();                 // velocity
+    std::vector<sdpmDouble> omega = _sol._pbl.getFrequencies(); // circular frequencies
+    // Geometry values
+    sdpmVector3d xRef = _sol._pbl.getRefCtr(); // reference center of rotation
+    // Organize panel and node data
+    sdpmMatrixXd cg = sdpmMatrixXd::Zero(_sol._nP, 3); // panel CG
+    sdpmMatrixXd s = sdpmVectorXd::Zero(_sol._nP);     // panel area
+    sdpmMatrixXd n = sdpmMatrixXd::Zero(_sol._nP, 3);  // panel compressible normal
+    sdpmMatrixXd nz = sdpmVectorXd::Zero(_sol._nP);    // panel z-normal
+    sdpmMatrixXd u = sdpmMatrixXd::Zero(_sol._nP, 3);  // steady velocity
+    sdpmVectorXd ux = sdpmVectorXd::Zero(_sol._nP);    // steady perturbation x-velocity
+    for (auto body : _sol._pbl.getBodies())
+    {
+        for (auto e : body->getElements())
+        {
+            cg.row(_sol._rows.at(e)) = e->getCg();
+            s(_sol._rows.at(e)) = e->getSurfaceArea();
+            n.row(_sol._rows.at(e)) = e->getCompressibleNormal();
+            nz(_sol._rows.at(e)) = e->getNormal()(2);
+            u.row(_sol._rows.at(e)) = _sol._u[e->getId() - 1];
+            ux(_sol._rows.at(e)) = _sol._u[e->getId() - 1](0) - uinf(0);
+        }
+    }
+
+    // Compute derivatives
+    std::cout << "Computing pressure derivatives... " << std::flush;
+    for (size_t ifq = 0; ifq < omega.size(); ++ifq)
+    {
+        // K matrices
+        sdpmMatrixXcd K = _sol._A[ifq].partialPivLu().inverse() * _sol._B[ifq];
+        std::vector<sdpmMatrixXcd> Kx(3);
+        for (size_t idim = 0; idim < 3; ++idim)
+            Kx[idim] = dN[idim] * L * K;
+
+        // Pressure stiffness terms (multiplied by 1)
+        // w.r.t angle of attack
+        sdpmVectorXcd cp_a = (2 / beta * u.col(0).array() - 2 * mach * mach / beta * ux.array()) * (uinf(0) * n.col(2).array() * n.col(0).array() + (Kx[0] * uinf(0) * n.col(2)).array()) + 2 * u.col(1).array() * (uinf(0) * n.col(2).array() * n.col(1).array() + (Kx[1] * uinf(0) * n.col(2)).array()) + 2 * u.col(2).array() * (-uinf(0) + uinf(0) * n.col(2).array() * n.col(2).array() + (Kx[2] * uinf(0) * n.col(2)).array());
+        _dCp1Re_a[ifq] = cp_a.real();
+        _dCp1Im_a[ifq] = cp_a.imag();
+        // w.r.t modal x-rotations
+        sdpmMatrixXcd cp_rx = (2 / beta * u.col(0).asDiagonal() - 2 * mach * mach / beta * ux.asDiagonal()) * (sdpmMatrixXd((uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) * n.col(0).asDiagonal() + Kx[0] * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) + 2 * u.col(1).asDiagonal() * (-uinf(2) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd((uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) * n.col(1).asDiagonal() + Kx[1] * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) + 2 * u.col(2).asDiagonal() * (uinf(1) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd((uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) * n.col(2).asDiagonal() + Kx[2] * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal());
+        _dCp1Re_rx[ifq] = cp_rx.real();
+        _dCp1Im_rx[ifq] = cp_rx.imag();
+        // w.r.t modal y-rotations
+        sdpmMatrixXcd cp_ry = -2 * u.col(0).asDiagonal() * uinf(2) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + (2 / beta * u.col(0).asDiagonal() - 2 * mach * mach / beta * ux.asDiagonal()) * (sdpmMatrixXd((uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) * n.col(0).asDiagonal() + Kx[0] * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) + 2 * u.col(1).asDiagonal() * (sdpmMatrixXd((uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) * n.col(1).asDiagonal() + Kx[1] * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) + 2 * u.col(2).asDiagonal() * (uinf(0) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd((uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) * n.col(2).asDiagonal() + Kx[2] * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal());
+        _dCp1Re_ry[ifq] = cp_ry.real();
+        _dCp1Im_ry[ifq] = cp_ry.imag();
+
+        // Pressure damping terms (multiplied by i*omega)
+        // w.r.t time derivative of angle of attack
+        sdpmVectorXcd cp_aDot = 2 * u.col(0).array() * (cg.col(2).array() - xRef(2)) + (2 / beta * u.col(0).array() - 2 * mach * mach / beta * ux.array()) * ((-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()) * n.col(0).array() + (Kx[0] * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array()) + 2 * u.col(1).array() * ((-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()) * n.col(1).array() + (Kx[1] * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array()) + 2 * u.col(2).array() * (-(cg.col(0).array() - xRef(0)) + (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()) * n.col(2).array() + (Kx[2] * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array()) + 2 * (K * uinf(0) * n.col(2)).array() * (1 - mach * mach * ux.array());
+        cp_aDot *= sdpmComplex(0., omega[ifq]);
+        _dCp1Re_aDot[ifq] = cp_aDot.real();
+        _dCp1Im_aDot[ifq] = cp_aDot.imag();
+        // w.r.t time derivative of z-displacement
+        sdpmVectorXcd cp_hDot = (2 / beta * u.col(0).array() - 2 * mach * mach / beta * ux.array()) * (n.col(2).array() * n.col(0).array() + (Kx[0] * n.col(2)).array()) + 2 * u.col(1).array() * (n.col(2).array() * n.col(1).array() + (Kx[1] * n.col(2)).array()) + 2 * u.col(2).array() * (-1 + n.col(2).array() * n.col(2).array() + (Kx[2] * n.col(2)).array());
+        cp_hDot *= sdpmComplex(0., omega[ifq]);
+        _dCp1Re_hDot[ifq] = cp_hDot.real();
+        _dCp1Im_hDot[ifq] = cp_hDot.imag();
+        // w.r.t time derivative of modal x-rotations
+        sdpmMatrixXcd cp_rxDot = sdpmComplex(0., 2 * omega[ifq]) * (sdpmMatrixXd::Identity(_sol._nP, _sol._nP) - sdpmMatrixXd(mach * mach * ux.asDiagonal())) * (K * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal());
+        _dCp1Re_rxDot[ifq] = cp_rxDot.real();
+        _dCp1Im_rxDot[ifq] = cp_rxDot.imag();
+        // w.r.t time derivative of modal y-rotations
+        sdpmMatrixXcd cp_ryDot = sdpmComplex(0., 2 * omega[ifq]) * (sdpmMatrixXd::Identity(_sol._nP, _sol._nP) - sdpmMatrixXd(mach * mach * ux.asDiagonal())) * (K * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal());
+        _dCp1Re_ryDot[ifq] = cp_ryDot.real();
+        _dCp1Im_ryDot[ifq] = cp_ryDot.imag();
+        // w.r.t time derivative of modal z-displacements
+        sdpmMatrixXcd cp_dzDot = (2 / beta * u.col(0).asDiagonal() - 2 * mach * mach / beta * ux.asDiagonal()) * (sdpmMatrixXd(n.col(2).asDiagonal()) * n.col(0).asDiagonal() + Kx[0] * n.col(2).asDiagonal()) + 2 * u.col(1).asDiagonal() * (sdpmMatrixXd(n.col(2).asDiagonal()) * n.col(1).asDiagonal() + Kx[1] * n.col(2).asDiagonal()) + 2 * u.col(2).asDiagonal() * (-sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd(n.col(2).asDiagonal()) * n.col(2).asDiagonal() + Kx[2] * n.col(2).asDiagonal());
+        cp_dzDot *= sdpmComplex(0., omega[ifq]);
+        _dCp1Re_dzDot[ifq] = cp_dzDot.real();
+        _dCp1Im_dzDot[ifq] = cp_dzDot.imag();
+
+        // Pressure inertia terms (multiplied by (i*omega)^2)
+        // w.r.t second time derivative of angle of attack
+        sdpmVectorXcd cp_aDotDot = -2 * omega[ifq] * omega[ifq] * (1 - mach * mach * ux.array()) * (K * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array();
+        _dCp1Re_aDotDot[ifq] = cp_aDotDot.real();
+        _dCp1Im_aDotDot[ifq] = cp_aDotDot.imag();
+        // w.r.t second time derivative of z-displacement
+        sdpmVectorXcd cp_hDotDot = -2 * omega[ifq] * omega[ifq] * (1 - mach * mach * ux.array()) * (K * n.col(2)).array();
+        _dCp1Re_hDotDot[ifq] = cp_hDotDot.real();
+        _dCp1Im_hDotDot[ifq] = cp_hDotDot.imag();
+        // w.r.t second time derivative of modal z-displacements
+        sdpmMatrixXcd cp_dzDotDot = -2 * omega[ifq] * omega[ifq] * (sdpmMatrixXd::Identity(_sol._nP, _sol._nP) - sdpmMatrixXd(mach * mach * ux.asDiagonal())) * (K * n.col(2).asDiagonal());
+        _dCp1Re_dzDotDot[ifq] = cp_dzDotDot.real();
+        _dCp1Im_dzDotDot[ifq] = cp_dzDotDot.imag();
+
+        // Load terms
+        sdpmMatrixXcd clz_rx = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_rx + cp_rxDot) * L1;
+        _dClz1Re_rx[ifq] = clz_rx.real();
+        _dClz1Im_rx[ifq] = clz_rx.imag();
+        sdpmMatrixXcd clz_ry = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_ry + cp_ryDot) * L1;
+        _dClz1Re_ry[ifq] = clz_ry.real();
+        _dClz1Im_ry[ifq] = clz_ry.imag();
+        sdpmMatrixXcd clz_dz = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_dzDot + cp_dzDotDot) * L1;
+        _dClz1Re_dz[ifq] = clz_dz.real();
+        _dClz1Im_dz[ifq] = clz_dz.imag();
+    }
+    std::cout << "done." << std::endl;
+}
+
+/**
+ * @brief Compute partial derivative of real part of GAF matrix entry at given frequency with respect to given mode
+ */
+std::vector<sdpmDouble> Gradient::computePartialsGafRe(size_t ifq, int imd, int irw, int jcl, std::string const &wrt)
+{
+    // Check input
+    sdpmMatrixXd dClz1;
+    if (wrt == "rx")
+        dClz1 = _dClz1Re_rx[ifq];
+    else if (wrt == "ry")
+        dClz1 = _dClz1Re_ry[ifq];
+    else if (wrt == "dz")
+        dClz1 = _dClz1Re_dz[ifq];
+    else
+        throw std::runtime_error("sdpm::Gradient::computePartialsGafRe: wrt parameter should be either 'rx', 'ry' or 'dz'!\n");
+
+    // Compute derivative
+    sdpmVectorXd dq = sdpmVectorXd::Zero(_nN);
+    if (jcl == imd)
+        dq += dClz1.transpose() * _dz[irw];
+    if (irw == imd && wrt == "dz")
+        dq -= _clz1Re[jcl][ifq];
+    return std::vector<sdpmDouble>(dq.data(), dq.data() + dq.size());
+}
+
+/**
+ * @brief Compute partial derivative of imaginary part of GAF matrix entry at given frequency with respect to given mode
+ */
+std::vector<sdpmDouble> Gradient::computePartialsGafIm(size_t ifq, int imd, int irw, int jcl, std::string const &wrt)
+{
+    // Check input
+    sdpmMatrixXd dClz1;
+    if (wrt == "rx")
+        dClz1 = _dClz1Im_rx[ifq];
+    else if (wrt == "ry")
+        dClz1 = _dClz1Im_ry[ifq];
+    else if (wrt == "dz")
+        dClz1 = _dClz1Im_dz[ifq];
+    else
+        throw std::runtime_error("sdpm::Gradient::computePartialsGafIm: wrt parameter should be either 'rx', 'ry' or 'dz'!\n");
+
+    // Compute derivative
+    sdpmVectorXd dq = sdpmVectorXd::Zero(_nN);
+    if (jcl == imd)
+        dq += dClz1.transpose() * _dz[irw];
+    if (irw == imd && wrt == "dz")
+        dq -= _clz1Im[jcl][ifq];
+    return std::vector<sdpmDouble>(dq.data(), dq.data() + dq.size());
+}
+
+/**
+ * @brief Test pressure derivatives by mulitplying with DOF and comparing to pressure
+ * @todo remove
+ */
+void Gradient::test(double a_amp, double h_amp)
+{
+    sdpmVectorXd cp1re = (_dCp1Re_a[0] + _dCp1Re_aDot[0] + _dCp1Re_aDotDot[0]) * a_amp + (_dCp1Re_hDot[0] + _dCp1Re_hDotDot[0]) * h_amp;
+    sdpmVectorXd cp1im = (_dCp1Im_a[0] + _dCp1Im_aDot[0] + _dCp1Im_aDotDot[0]) * a_amp + (_dCp1Im_hDot[0] + _dCp1Im_hDotDot[0]) * h_amp;
+    for (auto b : _sol._pbl.getBodies())
+    {
+        for (auto e : b->getElements())
+        {
+            sdpmDouble deltaCp1Re = _sol._cp1Re[0][0][e->getId() - 1] - cp1re(_sol._rows.at(e));
+            sdpmDouble deltaCp1Im = _sol._cp1Im[0][0][e->getId() - 1] - cp1im(_sol._rows.at(e));
+            if (abs(deltaCp1Re) > 1e-12 || abs(deltaCp1Im) > 1e-12)
+                throw std::runtime_error("sdpm::Gradient::test() failed!\n");
+        }
+    }
+}
+
+void Gradient::write(std::ostream &out) const
+{
+    out << "sdpm::Gradient on\n"
+        << _sol;
+}

sdpm/src/sdpmGradient.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 SDPMGRADIENT_H
+#define SDPMGRADIENT_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+#include <vector>
+#include <map>
+
+namespace sdpm
+{
+
+/**
+ * @brief Source-doublet panel gradient interface
+ * @authors Adrien Crovato
+ * @todo remove dummy test class once rigid motions (a and h) have been included in GAF matrix computation
+ */
+class SDPM_API Gradient : public sdpmObject
+{
+    Solver &_sol;                 ///< SDPM solver
+    size_t _nF;                   ///< number of frequencies
+    size_t _nM;                   ///< number of modes
+    size_t _nN;                   ///< number of nodes
+    std::map<Node *, int> _rowsN; ///< node row identification map
+
+    // Displacements and loads
+    std::vector<sdpmVectorXd> _dz;                  ///< z-displacement at nodes
+    std::vector<std::vector<sdpmVectorXd>> _clz1Re; ///< z-component of real part of pressure load at nodes
+    std::vector<std::vector<sdpmVectorXd>> _clz1Im; ///< z-component of real imaginary of pressure load at nodes
+
+    // Pressure derivatives (stiffness terms)
+    std::vector<sdpmVectorXd> _dCp1Re_a;  ///< gradient of real part of unsteady pressure w.r.t. angle of attack
+    std::vector<sdpmVectorXd> _dCp1Im_a;  ///< gradient of imaginary part of unsteady pressure w.r.t. angle of attack
+    std::vector<sdpmMatrixXd> _dCp1Re_rx; ///< gradient of real part of unsteady pressure w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Im_rx; ///< gradient of imaginary part of unsteady pressure w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Re_ry; ///< gradient of real part of unsteady pressure w.r.t. modal rotations about y
+    std::vector<sdpmMatrixXd> _dCp1Im_ry; ///< gradient of imaginary part of unsteady pressure w.r.t. modal rotations about y
+    // Pressure derivatives (damping terms)
+    std::vector<sdpmVectorXd> _dCp1Re_aDot;  ///< gradient of real part of unsteady pressure w.r.t. time derivative of angle of attack
+    std::vector<sdpmVectorXd> _dCp1Im_aDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of angle of attack
+    std::vector<sdpmVectorXd> _dCp1Re_hDot;  ///< gradient of real part of unsteady pressure w.r.t. time derivative of vertical displacement
+    std::vector<sdpmVectorXd> _dCp1Im_hDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of vertical displacement
+    std::vector<sdpmMatrixXd> _dCp1Re_rxDot; ///< gradient of real part of unsteady pressure w.r.t. time derivative of modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Im_rxDot; ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Re_ryDot; ///< gradient of real part of unsteady pressure w.r.t. time derivative of modal rotations about y
+    std::vector<sdpmMatrixXd> _dCp1Im_ryDot; ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of modal rotations about y
+    std::vector<sdpmMatrixXd> _dCp1Re_dzDot; ///< gradient of real part of unsteady pressure w.r.t. time derivative of modal displacements along z
+    std::vector<sdpmMatrixXd> _dCp1Im_dzDot; ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of modal displacements along z
+    // Pressure derivatives (inertia terms)
+    std::vector<sdpmVectorXd> _dCp1Re_aDotDot;  ///< gradient of real part of unsteady pressure w.r.t. second time derivative of angle of attack
+    std::vector<sdpmVectorXd> _dCp1Im_aDotDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. second time derivative of angle of attack
+    std::vector<sdpmVectorXd> _dCp1Re_hDotDot;  ///< gradient of real part of unsteady pressure w.r.t. second time derivative of vertical displacement
+    std::vector<sdpmVectorXd> _dCp1Im_hDotDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. second time derivative of vertical displacement
+    std::vector<sdpmMatrixXd> _dCp1Re_dzDotDot; ///< gradient of real part of unsteady pressure w.r.t. second time derivative of modal displacements along z
+    std::vector<sdpmMatrixXd> _dCp1Im_dzDotDot; ///< gradient of imaginary part of unsteady pressure w.r.t. second time derivative of modal displacements along z
+
+    // Loads derivatives
+    std::vector<sdpmMatrixXd> _dClz1Re_rx; ///< gradient of real part of unsteady loads z-component w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dClz1Im_rx; ///< gradient of imaginary part of unsteady loads z-component w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dClz1Re_ry; ///< gradient of real part of unsteady loads z-component w.r.t. modal rotations about y
+    std::vector<sdpmMatrixXd> _dClz1Im_ry; ///< gradient of imaginary part of unsteady loads z-component w.r.t. modal rotations about y
+    std::vector<sdpmMatrixXd> _dClz1Re_dz; ///< gradient of real part of unsteady loads z-component w.r.t. modal displacements along z
+    std::vector<sdpmMatrixXd> _dClz1Im_dz; ///< gradient of imaginary part of unsteady loads z-component w.r.t. modal displacements along z
+
+public:
+    Gradient(Solver &sol);
+    ~Gradient();
+
+    void run();
+    std::vector<sdpmDouble> computePartialsGafRe(size_t ifq, int imd, int irw, int jcl, std::string const &wrt);
+    std::vector<sdpmDouble> computePartialsGafIm(size_t ifq, int imd, int irw, int jcl, std::string const &wrt);
+
+    void test(double a_amp, double h_amp);
+
+#ifndef SWIG
+    virtual void write(std::ostream &out) const override;
+#endif
+};
+
+} // namespace sdpm
+
+#endif // SDPMGRADIENT_H

sdpm/src/sdpmMotion.cpp

@@ -23,7 +23,7 @@ using namespace sdpm;
 Motion::Motion(size_t n) : sdpmObject(),
                            _aAmp(0.), _hAmp(0.),
                            _xRef(0.), _zRef(0.),
-                           _mAmp(0.), _dzMod(n, 0.), _rxMod(n, 0.), _ryMod(n, 0.)
+                           _dzMod(n, 0.), _rxMod(n, 0.), _ryMod(n, 0.)
 {
 }
 
@@ -43,12 +43,11 @@ void Motion::setRigid(double aAmp, double hAmp, double xRef, double zRef)
  */
 void Motion::setFlexible(double mAmp, std::vector<std::vector<double>> const &xMod, std::vector<Node *> const &nodes)
 {
-    _mAmp = mAmp;
     for (size_t i = 0; i < nodes.size(); ++i)
     {
-        _dzMod[nodes[i]->getId() - 1] = xMod[i][0];
-        _rxMod[nodes[i]->getId() - 1] = xMod[i][1];
-        _ryMod[nodes[i]->getId() - 1] = xMod[i][2];
+        _dzMod[nodes[i]->getId() - 1] = mAmp * xMod[i][0];
+        _rxMod[nodes[i]->getId() - 1] = mAmp * xMod[i][1];
+        _ryMod[nodes[i]->getId() - 1] = mAmp * xMod[i][2];
     }
 }
 
@@ -67,7 +66,7 @@ sdpmVector3d Motion::computeSteady(Element const &e, sdpmVector3d const &ui)
     {
         sdpmDouble rx = _rxMod[n->getId() - 1];
         sdpmDouble ry = _ryMod[n->getId() - 1];
-        flex += _mAmp * sdpmVector3d(ui(2) * ry, ui(2) * rx, -ui(0) * ry - ui(1) * rx) / e.getNodes().size();
+        flex += sdpmVector3d(ui(2) * ry, ui(2) * rx, -ui(0) * ry - ui(1) * rx) / e.getNodes().size();
     }
     return rigd + flex;
 }
@@ -78,11 +77,11 @@ sdpmVector3d Motion::computeSteady(Element const &e, sdpmVector3d const &ui)
 sdpmVector3d Motion::computeHarmonic(Element const &e)
 {
     // Rigid contribution
-    sdpmVector3d cg = e.getCompressibleCg();
+    sdpmVector3d cg = e.getCg();
     sdpmVector3d rigd(-_aAmp * (cg(2) - _zRef), 0., _aAmp * (cg(0) - _xRef) + _hAmp);
     // Flexible contribution
     sdpmVector3d flex = sdpmVector3d::Zero();
     for (auto n : e.getNodes())
-        flex += _mAmp * sdpmVector3d(0., 0., _dzMod[n->getId() - 1]) / e.getNodes().size();
+        flex += sdpmVector3d(0., 0., _dzMod[n->getId() - 1]) / e.getNodes().size();
     return rigd + flex;
 }

sdpm/src/sdpmMotion.h

@@ -34,14 +34,13 @@ class SDPM_API Motion : public sdpmObject
     sdpmDouble _hAmp;               ///< plunge amplitude
     sdpmDouble _xRef;               ///< x-coordinate of center of pitch
     sdpmDouble _zRef;               ///< z-coordinate of center of pitch
-    sdpmDouble _mAmp;               ///< modal amplitude
     std::vector<sdpmDouble> _dzMod; ///< modal displacements along z
     std::vector<sdpmDouble> _rxMod; ///< modal rotations about x
     std::vector<sdpmDouble> _ryMod; ///< modal rotations about y
 
 public:
     Motion(size_t n);
-    ~Motion(){};
+    ~Motion() {};
 
     sdpmDouble getPitchAmplitude() const { return _aAmp; }
     std::vector<sdpmDouble> const &getModeZ() const { return _dzMod; }

sdpm/src/sdpmProblem.h

@@ -31,7 +31,8 @@ namespace sdpm
  */
 class SDPM_API Problem : public sdpmObject
 {
-    friend class Adjoint; // so that Adjoint can access the problem inputs
+    friend class Adjoint;  // so that Adjoint can access the problem inputs
+    friend class Gradient; // so that Gradient can access the problem outputs
 
     // Mesh
     Mesh &_msh;                  ///< mesh data structure

sdpm/src/sdpmSolver.cpp

@@ -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

sdpm/src/sdpmSolver.h

@@ -31,7 +31,8 @@ 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

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., xf, 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., xf, 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 math
+
+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 * math.pi / 180
+    kred = 0.133 # reduced frequency
+    a_amp = 0.5 * 0.25 * math.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, xf, 0.)
+    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)
+
+    # test
+    grd.test(a_amp, h_amp) # TODO, replace by GAF matrix gradient check once implemented
+    tests = CTests()
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()

sdpm/tests/lann_xsonic.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., xf, 0.)
     dcp = interpolate_pressure(build_fpath('models/lann_dcp_rans.csv', __file__, 1), wing.getNodes())
     wing.setTransonicPressureGrad(dcp[:,0])
     pbl.addBody(wing)
@@ -78,14 +78,14 @@ def main():
     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.45, 5e-2))
-    tests.add(CTest('Ang(CL1)', cmath.phase(cl1) * 180 / math.pi, -3.08, 5e-2))
-    tests.add(CTest('Abs(CD1)', abs(cd1) / amp / math.pi, 0.193, 5e-2))
-    tests.add(CTest('Ang(CD1)', cmath.phase(cd1) * 180 / math.pi, 4.2, 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, -12.3, 5e-2))
-    tests.add(CTest('Abs(CM1)', abs(cm1) / amp / math.pi, 0.51, 5e-2))
-    tests.add(CTest('Ang(CM1)', cmath.phase(cm1) * 180 / math.pi, -172.3, 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