diff --git a/CMakeLists.txt b/CMakeLists.txt
index 105e3ac1a04e741e94011702086ba31f1df61ab1..82a5848331bdf4e0a5d841461bdc230c69593693 100644
--- a/CMakeLists.txt
+++ b/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
diff --git a/sdpm/_src/sdpmw.i b/sdpm/_src/sdpmw.i
index b27441a228f732ce52aa2a230dbc9b2b332019ee..289bf69eecf24636b12a54e54daffe7eafe065e2 100644
--- a/sdpm/_src/sdpmw.i
+++ b/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"
diff --git a/sdpm/src/sdpm.h b/sdpm/src/sdpm.h
index 0fa95980af2d93c94b558621de308d04a43c3003..70878202c7b043945413b18d8abc993ce878054e 100644
--- a/sdpm/src/sdpm.h
+++ b/sdpm/src/sdpm.h
@@ -81,6 +81,7 @@ class Builder;
class Solver;
class SolverTransonic;
class Adjoint;
+class Gradient;
} // namespace sdpm
#endif // SDPM_H
diff --git a/sdpm/src/sdpmGradient.cpp b/sdpm/src/sdpmGradient.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cfa0a806155ec5f09cbf626009ffe9eac5cca13f
--- /dev/null
+++ b/sdpm/src/sdpmGradient.cpp
@@ -0,0 +1,310 @@
+/*
+ * 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;
+}
diff --git a/sdpm/src/sdpmGradient.h b/sdpm/src/sdpmGradient.h
new file mode 100644
index 0000000000000000000000000000000000000000..4da66726d476d1d86ae45eb6305c7d94d195b8a6
--- /dev/null
+++ b/sdpm/src/sdpmGradient.h
@@ -0,0 +1,97 @@
+/*
+ * 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
diff --git a/sdpm/src/sdpmProblem.h b/sdpm/src/sdpmProblem.h
index 0b7022e37bf8f4f3fa71c45546d4abae460cbb61..046cfe41ab27746314c6883496c6b7636a8d4867 100644
--- a/sdpm/src/sdpmProblem.h
+++ b/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
diff --git a/sdpm/src/sdpmSolver.h b/sdpm/src/sdpmSolver.h
index ab5cc32eb10585b3df7a5e30664b7230fff81ba6..fa2c5ff833e8d0f5da3e8ad03741d31e146cbfb8 100644
--- a/sdpm/src/sdpmSolver.h
+++ b/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
diff --git a/sdpm/tests/agard_grad.py b/sdpm/tests/agard_grad.py
new file mode 100644
index 0000000000000000000000000000000000000000..7d06d44285d17337fe8a6167dea63c92684ab11b
--- /dev/null
+++ b/sdpm/tests/agard_grad.py
@@ -0,0 +1,103 @@
+# -*- 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()
diff --git a/sdpm/tests/lann_grad.py b/sdpm/tests/lann_grad.py
new file mode 100644
index 0000000000000000000000000000000000000000..85d5e335932f8104a71b333c9bda59c2e03c8b0a
--- /dev/null
+++ b/sdpm/tests/lann_grad.py
@@ -0,0 +1,83 @@
+# -*- 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()