USCSDPM v1.0

Version 1.0 - Unsteady compressible source and doublet panel method

USCSDPM, shortened SDPM, is an implementation of the original source and doublet panel method proposed by Morino (L. Morino and P. Chen, Steady and oscillatory subsonic and supersonic aerodynamics around complex configurations, AIAA Journal, 1975) and revisited by Dimitriadis (M. Martínez and G. Dimitriadis, Subsonic source and doublet panel methods, Journal of Fluids and Structures, 2022). The main purpose of SDPM is to provide steady and unsteady pressure distributions on fixed and moving wings. The code is currently capable of computing steady or unsteady, compressible or incompressible flows on isolated or several wing surfaces.

Features

• C++ code interfaced in Python using SWIG
• Steady or unsteady flows caused by forced harmonic pitching motion of the boundary
• Incompressible or compressible flows
• Isolated or multiple wings

Limitations

• Flow are currently limited to the compressible subsonic, but not transonic regime
• No gradient computation method is currently featured
• Only rigid pitching motion can be currently imposed on the moving boundary

These limitations will be addressed in future developments.

Related issues

• #4 (closed) implemented
• #5 (closed) fixed

Tests

• Ubuntu 22.04 (gcc 11, python 3.10)
• Windows 10 (vs 2022, python 3.8)

sdpm/src/sdpmSolver.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 "sdpmSolver.h"
+#include "sdpmBuilder.h"
+#include "sdpmProblem.h"
+#include "sdpmBody.h"
+#include "sdpmMesh.h"
+#include "sdpmGmshExport.h"
+#include "sdpmResults.h"
+#include "sdpm_extras.h"
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+
+using namespace sdpm;
+
+Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
+{
+    // Say hi
+    std::cout << std::endl;
+    std::cout << "*******************************************************************************" << std::endl;
+    std::cout << "** Hi! My name is SDPM v1.0-2301                                             **" << std::endl;
+    std::cout << "** ULiege 2023-2023                                                          **" << std::endl;
+    std::cout << "*******************************************************************************" << std::endl;
+
+    // Check the problem
+    _pbl.check();
+
+    // Count panel and set up global->local id map
+    _nP = 0;
+    int i = 0;
+    for (auto body : _pbl.getBodies())
+    {
+        _nP += body->getElements().size();
+        for (auto e : body->getElements())
+            _rows[e] = i++;
+    }
+
+    // Get sizes
+    size_t nm = _pbl.getNumModes();
+    size_t nf = _pbl.getFrequencies().size();
+    size_t ne = _pbl.getElements().size();
+
+    // Resize AIC matrics
+    _A0 = sdpmMatrixXd(_nP, _nP);
+    _B0 = sdpmMatrixXd(_nP, _nP);
+    _A.resize(nf, sdpmMatrixXcd::Zero(_nP, _nP));
+    _B.resize(nf, sdpmMatrixXcd::Zero(_nP, _nP));
+
+    // Set up variables
+    _phi.resize(ne, 0.);
+    _tau.resize(ne, 0.);
+    _mu.resize(ne, 0.);
+    _u.resize(ne, sdpmVector3d::Zero());
+    _cp.resize(ne, 0.);
+    _u1.resize(nm, std::vector<std::vector<sdpmVector3cd>>(nf, std::vector<sdpmVector3cd>(ne, sdpmVector3cd::Zero())));
+    _cp1.resize(nm, std::vector<std::vector<sdpmComplex>>(nf, std::vector<sdpmComplex>(ne, sdpmComplex(0., 0.))));
+    _cl1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+    _cd1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+    _cs1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+    _cm1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+
+    // Display configuration
+    std::cout << "--- Physical groups ---\n";
+    std::cout << "Number of panels: " << _nP << "\n";
+    for (auto b : _pbl.getBodies())
+        std::cout << *b << "\n";
+    std::cout << "--- Freestream conditions ---\n"
+              << std::setprecision(5)
+              << "Angle of attack: " << _pbl.getAoa() * 180 / 3.14159 << " deg\n"
+              << "Angle of sideslip: " << _pbl.getAos() * 180 / 3.14159 << " deg\n"
+              << "Mach number: " << _pbl.getMach() << "\n"
+              << "Velocity: " << _pbl.getVelocity().norm() << "\n";
+    std::cout << "--- Reference data ---\n"
+              << "Surface area: " << _pbl.getRefArea() << "\n"
+              << "Chord length: " << _pbl.getRefLgth() << "\n"
+              << "Origin: (" << _pbl.getRefCtr().transpose() << ")\n";
+    std::cout << "--- Unsteady parameters ---\n"
+              << "Number of frequencies: " << _pbl.getFrequencies().size() << "\n"
+              << "Number of motions: " << _pbl.getNumModes() << "\n"
+              << std::endl;
+}
+
+Solver::~Solver()
+{
+    std::cout << "~sdpm::Solver()\n";
+}
+
+/**
+ * @brief Update elements and AIC matrices
+ */
+void Solver::update()
+{
+    // Update elements and get values
+    _pbl.updateElements();
+    std::vector<double> omega = _pbl.getFrequencies();
+    sdpmDouble ui = _pbl.getVelocity().norm();
+    sdpmDouble mi = _pbl.getMach();
+    sdpmDouble beta = _pbl.getCompressibility();
+    sdpmComplex img(0, 1);
+
+    // Compute AIC
+    std::cout << "Computing AIC matrices... " << std::flush;
+    std::vector<Body *> bodies = _pbl.getBodies();
+    _A0.setZero();
+    _B0.setZero();
+    std::fill(_A.begin(), _A.end(), sdpmMatrixXcd::Zero(_nP, _nP));
+    std::fill(_B.begin(), _B.end(), sdpmMatrixXcd::Zero(_nP, _nP));
+    for (auto jbody : bodies)
+    {
+        // body
+        for (auto ej : jbody->getElements())
+        {
+            sdpmDouble nx = ej->getCompressibleNormal()(0);                /// x-component of normal vector
+            sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ej);      // global to local transformation matrix
+            sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(false, ej); // coordinates of source panel in local axes
+            for (auto ibody : bodies)
+            {
+                for (auto ei : ibody->getElements())
+                {
+                    std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ej, sCoords, g2l); // steady AIC
+                    _A0(_rows.at(ei), _rows.at(ej)) = mutau[0];
+                    _B0(_rows.at(ei), _rows.at(ej)) = mutau[1];
+                    for (size_t fk = 0; fk < omega.size(); ++fk)
+                    {
+                        sdpmDouble om = omega[fk] * mi / ui / beta;
+                        sdpmVector3d r = Builder::computeDistVec(false, ei, ej);
+                        sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
+                        _A[fk](_rows.at(ei), _rows.at(ej)) = (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mi * nx) * exp;
+                        _B[fk](_rows.at(ei), _rows.at(ej)) = mutau[1] * exp;
+                    }
+                }
+            }
+        }
+        // wake
+        Wake const &wake = jbody->getWake();
+        std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getElMap();
+        std::map<Element *, double> lagMap = wake.getLagMap();
+        for (auto ew : wake.getElements())
+        {
+            sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ew);      // global to local transformation matrix
+            sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(false, ew); // coordinates of source panel in local axes
+            for (auto ibody : bodies)
+            {
+                for (auto ei : ibody->getElements())
+                {
+                    std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ew, sCoords, g2l); // steady AIC
+                    _A0(_rows.at(ei), _rows.at(wwMap.at(ew).first)) += mutau[0];
+                    _A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= mutau[0];
+                    for (size_t fk = 0; fk < omega.size(); ++fk)
+                    {
+                        sdpmDouble om = omega[fk] * mi / ui / beta;
+                        sdpmVector3d r = Builder::computeDistVec(false, ei, ew);
+                        sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
+                        sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
+                        _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) += waktrm;
+                        _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= waktrm;
+                    }
+                }
+            }
+        }
+    }
+    // symmetry
+    if (_pbl.getSym())
+    {
+        for (auto jbody : bodies)
+        {
+            // body
+            for (auto ej : jbody->getElements())
+            {
+                sdpmDouble nx = ej->getCompressibleNormal()(0);               /// x-component of normal vector
+                sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ej);      // global to local transformation matrix
+                sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(true, ej); // coordinates of source panel in local axes
+                for (auto ibody : bodies)
+                {
+                    for (auto ei : ibody->getElements())
+                    {
+                        std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ej, sCoords, g2l); // steady AIC
+                        _A0(_rows.at(ei), _rows.at(ej)) -= mutau[0];
+                        _B0(_rows.at(ei), _rows.at(ej)) -= mutau[1];
+                        for (size_t fk = 0; fk < omega.size(); ++fk)
+                        {
+                            sdpmDouble om = omega[fk] * mi / ui / beta;
+                            sdpmVector3d r = Builder::computeDistVec(true, ei, ej);
+                            sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
+                            _A[fk](_rows.at(ei), _rows.at(ej)) -= (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mi * nx) * exp;
+                            _B[fk](_rows.at(ei), _rows.at(ej)) -= mutau[1] * exp;
+                        }
+                    }
+                }
+            }
+            // wake
+            Wake const &wake = jbody->getWake();
+            std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getElMap();
+            std::map<Element *, double> lagMap = wake.getLagMap();
+            for (auto ew : wake.getElements())
+            {
+                sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ew);      // global to local transformation matrix
+                sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(true, ew); // coordinates of source panel in local axes
+                for (auto ibody : bodies)
+                {
+                    for (auto ei : ibody->getElements())
+                    {
+                        std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ew, sCoords, g2l); // steady AIC
+                        _A0(_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= mutau[0];
+                        _A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) += mutau[0];
+                        for (size_t fk = 0; fk < omega.size(); ++fk)
+                        {
+                            sdpmDouble om = omega[fk] * mi / ui / beta;
+                            sdpmVector3d r = Builder::computeDistVec(true, ei, ew);
+                            sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
+                            sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
+                            _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= waktrm;
+                            _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) += waktrm;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    std::cout << "done." << std::endl;
+}
+
+/**
+ * @brief Solve the steady potential equation
+ */
+void Solver::run()
+{
+    std::cout << "Running steady SDPM solver... " << std::endl;
+    std::cout << "computing sources... " << std::flush;
+    sdpmVector3d ui = _pbl.getVelocity(); // freestream velocity
+    sdpmVectorXd etau(_nP);
+    for (auto body : _pbl.getBodies())
+        for (auto e : body->getElements())
+            etau(_rows[e]) = ui.dot(e->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
+    std::cout << "done." << std::endl;
+
+    std::cout << "solving for doublets... " << std::flush;
+    sdpmVectorXd emu(_nP);
+    emu = Eigen::PartialPivLU<sdpmMatrixXd>(_A0).solve(_B0 * etau);
+    std::cout << "done." << std::endl;
+
+    std::cout << "computing flow and loads on bodies... " << std::flush;
+    post(etau, emu);
+    std::cout << "done." << std::endl
+              << std::endl;
+}
+
+/**
+ * @brief Solve the unsteady potential equation
+ * @todo remove output
+ */
+void Solver::runUnsteady()
+{
+    std::cout << "Running unsteady SDPM solver... " << std::endl;
+    for (size_t imd = 0; imd < _pbl.getNumModes(); ++imd)
+    {
+        // Compute motion-induced velocity
+        for (auto body : _pbl.getBodies())
+            body->computeVelocity(imd);
+        // Solve
+        for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
+        {
+            std::cout << "Mode #" << imd << ", frequency #" << ifq << std::endl;
+            std::cout << "computing sources... " << std::flush;
+            sdpmVectorXcd etau(_nP);
+            for (auto body : _pbl.getBodies())
+            {
+                std::vector<sdpmVector3cd> uc0 = body->getSteadyVelocity(imd);
+                std::vector<sdpmVector3cd> uc1 = body->getHarmonicVelocity(imd);
+                std::vector<Element *> elems = body->getElements();
+                for (size_t i = 0; i < elems.size(); ++i)
+                    etau(_rows[elems[i]]) = (uc0[i] + uc1[i] * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]).conjugate().dot(elems[i]->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
+            }
+            std::cout << "done." << std::endl;
+
+            std::cout << "solving for doublets... " << std::flush;
+            sdpmVectorXcd emu(_nP);
+            emu = Eigen::PartialPivLU<sdpmMatrixXcd>(_A[ifq]).solve(_B[ifq] * etau);
+            std::cout << "done." << std::endl;
+
+            std::cout << "computing flow and loads on bodies... " << std::flush;
+            postUnsteady(imd, ifq, etau, emu);
+            std::cout << "done." << std::endl;
+        }
+    }
+    std::cout << std::endl;
+}
+
+/**
+ * @brief Write the results
+ */
+void Solver::save(GmshExport const &writer, std::string const &suffix)
+{
+    // Get data
+    size_t nE = _pbl.getElements().size();
+    size_t nM = _pbl.getNumModes();
+    std::vector<sdpmDouble> omega = _pbl.getFrequencies();
+    size_t nF = omega.size();
+    sdpmDouble toDeg = 180. / 3.14159;
+
+    // Write files
+    std::cout << "Saving files... " << std::endl;
+    // set up results
+    Results results;
+    // steady
+    results.scalarsAtElems["phi"] = &_phi;
+    results.scalarsAtElems["mu"] = &_mu;
+    results.scalarsAtElems["tau"] = &_tau;
+    results.vectorsAtElems["u"] = &_u;
+    results.scalarsAtElems["cp"] = &_cp;
+    // unsteady
+    std::vector<std::vector<std::vector<sdpmVector3d>>> u1_abs(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(nE)));
+    std::vector<std::vector<std::vector<sdpmVector3d>>> u1_arg(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(nE)));
+    std::vector<std::vector<std::vector<sdpmDouble>>> cp1_abs(nM, std::vector<std::vector<sdpmDouble>>(nF, std::vector<sdpmDouble>(nE)));
+    std::vector<std::vector<std::vector<sdpmDouble>>> cp1_arg(nM, std::vector<std::vector<sdpmDouble>>(nF, std::vector<sdpmDouble>(nE)));
+    for (size_t imd = 0; imd < nM; ++imd)
+    {
+        for (size_t ifq = 0; ifq < nF; ++ifq)
+        {
+            // set names
+            std::stringstream u1_abs_name, u1_arg_name, cp1_abs_name, cp1_arg_name;
+            u1_abs_name << "|u1|_" << imd << "_" << ifq;
+            u1_arg_name << "< u1_" << imd << "_" << ifq;
+            cp1_abs_name << "|cp1|_" << imd << "_" << ifq;
+            cp1_arg_name << "< cp1_" << imd << "_" << ifq;
+            // set values
+            for (size_t i = 0; i < nE; ++i)
+            {
+                u1_abs[imd][ifq][i] = _u1[imd][ifq][i].array().abs();
+                u1_arg[imd][ifq][i] = _u1[imd][ifq][i].array().arg() * toDeg;
+                cp1_abs[imd][ifq][i] = std::abs(_cp1[imd][ifq][i]);
+                cp1_arg[imd][ifq][i] = std::arg(_cp1[imd][ifq][i]) * toDeg;
+            }
+            results.vectorsAtElems[u1_abs_name.str()] = &u1_abs[imd][ifq];
+            results.vectorsAtElems[u1_arg_name.str()] = &u1_arg[imd][ifq];
+            results.scalarsAtElems[cp1_abs_name.str()] = &cp1_abs[imd][ifq];
+            results.scalarsAtElems[cp1_arg_name.str()] = &cp1_arg[imd][ifq];
+        }
+    }
+    // save results
+    writer.save(results, suffix);
+    // save aerodynamic loads breakdown
+    std::cout << "writing file: aeroloads_breakdown.dat... " << std::flush;
+    std::ofstream outfile;
+    outfile.open("aeroloads_breakdown.dat");
+    for (auto b : _pbl.getBodies())
+    {
+        // steady
+        outfile << "# Body - " << b->getName() << " (" << b->getElements().size() << " elements)" << std::endl;
+        outfile << std::fixed << std::setprecision(5)
+                << "Cl = " << b->getLiftCoeff() << "\n"
+                << "Cd = " << b->getDragCoeff() << "\n"
+                << "Cs = " << b->getSideforceCoeff() << "\n"
+                << "Cm = " << b->getMomentCoeff() << "\n";
+        // unsteady
+        for (size_t imd = 0; imd < nM; ++imd)
+        {
+            for (size_t ifq = 0; ifq < nF; ++ifq)
+            {
+                outfile << "|Cl1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadyLiftCoeff(imd, ifq)) << "\n"
+                        << "< Cl1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadyLiftCoeff(imd, ifq)) * toDeg << "\n"
+                        << "|Cd1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadyDragCoeff(imd, ifq)) << "\n"
+                        << "< Cd1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadyDragCoeff(imd, ifq)) * toDeg << "\n"
+                        << "|Cs1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadySideforceCoeff(imd, ifq)) << "\n"
+                        << "< Cs1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadySideforceCoeff(imd, ifq)) * toDeg << "\n"
+                        << "|Cm1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadyMomentCoeff(imd, ifq)) << "\n"
+                        << "< Cm1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadyMomentCoeff(imd, ifq)) * toDeg << "\n"
+                        << std::endl;
+            }
+        }
+    }
+    outfile.close();
+    std::cout << "done." << std::endl;
+}
+
+/**
+ * @brief Compute steady flow functionals and loads on bodies
+ */
+void Solver::post(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
+{
+    // Store source and doublet strength and compute perturbation potential on elements
+    sdpmVectorXd ephi = -_A0 * emu - _B0 * etau;
+    std::vector<Body *> bodies = _pbl.getBodies();
+    for (auto body : bodies)
+    {
+        for (auto e : body->getElements())
+        {
+            size_t id = e->getId() - 1;
+            int i = _rows[e];
+            _tau[id] = etau(i);
+            _mu[id] = emu(i);
+            _phi[id] = ephi(i);
+        }
+    }
+
+    // Transfer doublet at nodes
+    std::vector<sdpmDouble> muNodes(_pbl.getNodes().size(), 0);
+    for (auto b : bodies)
+    {
+        std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+        for (auto n : b->getNodes())
+            for (auto e : neMap.at(n))
+                muNodes[n->getId() - 1] += _mu[e->getId() - 1] / neMap.at(n).size();
+    }
+
+    // Compute surface flow functionals
+    sdpmVector3d ui = _pbl.getVelocity(); // freestream velocity
+    sdpmDouble mi = _pbl.getMach();       // freestream speed of sound
+    for (auto body : bodies)
+    {
+        for (auto e : body->getElements())
+        {
+            sdpmVector3d u = ui - (_tau[e->getId() - 1] * e->getCompressibleNormal() + e->computeGradient(muNodes)).cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.));
+            _u[e->getId() - 1] = u;                                                                        // velocity
+            _cp[e->getId() - 1] = 1 - (u.dot(u) - mi * mi * (u(0) - ui(0)) * (u(0) - ui(0))) / ui.dot(ui); // nonlinear pressure coefficient
+        }
+    }
+
+    // Compute aerodynamic loads
+    sdpmDouble sRef = _pbl.getRefArea();      // reference surface area
+    sdpmDouble cRef = _pbl.getRefLgth();      // reference length
+    sdpmVector3d xRef = _pbl.getRefCtr();     // reference center
+    sdpmVector3d dragDir = _pbl.getDragDir(); // drag direction
+    sdpmVector3d sideDir = _pbl.getSideDir(); // sideforce direction
+    sdpmVector3d liftDir = _pbl.getLiftDir(); // lift direction
+    _cl = 0;
+    _cd = 0;
+    _cs = 0;
+    _cm = 0;
+    for (auto b : _pbl.getBodies())
+    {
+        // Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element and average at node
+        std::vector<Node *> nods = b->getNodes();
+        std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+        std::vector<sdpmVector3d> nLoads(nods.size(), sdpmVector3d::Zero());
+        for (size_t i = 0; i < nods.size(); ++i)
+            for (auto e : neMap.at(nods[i]))
+                nLoads[i] += -_cp[e->getId() - 1] * e->getSurfaceArea() * e->getNormal() / e->getNodes().size();
+        b->setLoads(nLoads);
+
+        // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
+        // force coefficients along x, y and z directions and pitching moment around y axis
+        sdpmDouble bcm = 0;
+        sdpmVector3d cf(0, 0, 0);
+        for (size_t i = 0; i < nods.size(); ++i)
+        {
+            sdpmVector3d l = nods[i]->getCoords() - xRef; // lever arm
+            sdpmVector3d cfi = nLoads[i];                 // normalized force
+            cf += cfi;
+            bcm += (l(2) * cfi(0) - l(0) * cfi(2)) / (sRef * cRef); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+        }
+        b->setMomentCoeff(bcm);
+        // rotate to flow direction
+        b->setDragCoeff(cf.dot(dragDir) / sRef);
+        b->setSideforceCoeff(cf.dot(sideDir) / sRef);
+        b->setLiftCoeff(cf.dot(liftDir) / sRef);
+        // compute total
+        _cl += b->getLiftCoeff();
+        _cd += b->getDragCoeff();
+        _cs += b->getSideforceCoeff();
+        _cm += bcm;
+    }
+}
+
+/**
+ * @brief Compute unsteady flow functionals and loads on bodies
+ */
+void Solver::postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu)
+{
+    // Transfer doublet at nodes
+    std::vector<sdpmDouble> muNodesRe(_pbl.getNodes().size(), 0), muNodesIm(_pbl.getNodes().size(), 0);
+    for (auto b : _pbl.getBodies())
+    {
+        std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+        for (auto n : b->getNodes())
+        {
+            for (auto e : neMap.at(n))
+            {
+                muNodesRe[n->getId() - 1] += emu(_rows.at(e)).real() / neMap.at(n).size();
+                muNodesIm[n->getId() - 1] += emu(_rows.at(e)).imag() / neMap.at(n).size();
+            }
+        }
+    }
+
+    // Compute flow and loads
+    sdpmVector3d ui = _pbl.getVelocity();
+    sdpmDouble mi = _pbl.getMach();
+    sdpmDouble omega = _pbl.getFrequencies()[ifq];
+    sdpmDouble sRef = _pbl.getRefArea();
+    sdpmDouble cRef = _pbl.getRefLgth();
+    sdpmVector3d xRef = _pbl.getRefCtr();
+    sdpmDouble ca = cos(_pbl.getAos());
+    sdpmDouble sa = sin(_pbl.getAos());
+    _cl1[imd][ifq] = sdpmComplex(0., 0.);
+    _cd1[imd][ifq] = sdpmComplex(0., 0.);
+    _cs1[imd][ifq] = sdpmComplex(0., 0.);
+    _cm1[imd][ifq] = sdpmComplex(0., 0.);
+    for (auto b : _pbl.getBodies())
+    {
+        // compute velocity
+        std::vector<sdpmVector3cd> uc0 = b->getSteadyVelocity(imd);
+        std::vector<sdpmVector3cd> uc1 = b->getHarmonicVelocity(imd);
+        std::vector<Element *> elems = b->getElements();
+        for (size_t i = 0; i < elems.size(); ++i)
+            _u1[imd][ifq][elems[i]->getId() - 1] = uc0[i] + uc1[i] * sdpmComplex(0., 1.) * omega - (etau[_rows[elems[i]]] * elems[i]->getCompressibleNormal() + (elems[i]->computeGradient(muNodesRe) + sdpmComplex(0, 1) * elems[i]->computeGradient(muNodesIm))).cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.));
+        // compute pressure
+        std::vector<std::vector<sdpmComplex>> cp012(elems.size(), std::vector<sdpmComplex>(3));
+        for (size_t i = 0; i < elems.size(); ++i)
+        {
+            // phi_t, phi_x
+            sdpmVector3cd phit(std::conj(-sdpmComplex(0, 1) * omega * emu(_rows[elems[i]])), 0., -sdpmComplex(0, 1) * omega * emu(_rows[elems[i]]));
+            sdpmVector3cd phix(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(0), _u[elems[i]->getId() - 1](0) - ui(0), _u1[imd][ifq][elems[i]->getId() - 1](0));
+            // 1 - phi_t
+            sdpmVectorXcd ecp(5);
+            ecp << 0., -2. * phit(0) / ui.dot(ui), 1., -2. * phit(2) / ui.dot(ui), 0.;
+            // phi_xx, phi_tt, phi_tx
+            ecp += mi * mi / ui.dot(ui) * (convolve(phix, phix) + convolve(phit, phit) + convolve(phix, phit) / ui.norm());
+            // conv(u_k,u_k)
+            for (size_t j = 0; j < 3; ++j)
+            {
+                sdpmVector3cd utrm(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(j), _u[elems[i]->getId() - 1](j), _u1[imd][ifq][elems[i]->getId() - 1](j));
+                ecp -= convolve(utrm, utrm) / ui.dot(ui);
+            }
+            // extract 0th, 1st and 2nd harmonics
+            for (size_t u = 0; u < 3; ++u)
+                cp012[i][u] = ecp(2 + u);
+            _cp1[imd][ifq][elems[i]->getId() - 1] = cp012[i][1];
+        }
+
+        // compute loads on nodes
+        std::vector<Node *> nods = b->getNodes();
+        std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+        std::vector<std::vector<sdpmVector3cd>> nLoads012(nods.size(), std::vector<sdpmVector3cd>(3, sdpmVector3cd::Zero()));
+        for (size_t i = 0; i < nods.size(); ++i)
+            for (auto e : neMap.at(nods[i]))
+                for (size_t u = 0; u < 3; ++u)
+                    nLoads012[i][u] += -cp012[_rows.at(e)][u] * e->getSurfaceArea() * e->getNormal() / e->getNodes().size();
+        std::vector<sdpmVector3cd> nLoads1(nods.size());
+        for (size_t i = 0; i < nods.size(); ++i)
+            nLoads1[i] = nLoads012[i][1];
+        b->setUnsteadyLoads(imd, ifq, nLoads1);
+
+        // compute moment coefficient and total load on body
+        sdpmComplex bcm1(0., 0.);
+        std::vector<sdpmVector3cd> cf(3, sdpmVector3cd::Zero());
+        for (size_t i = 0; i < nods.size(); ++i)
+        {
+            sdpmVector3d l = nods[i]->getCoords() - xRef;
+            for (size_t u = 0; u < 3; ++u)
+                cf[u] += nLoads012[i][u];
+            bcm1 += (l(2) * nLoads012[i][1](0) - l(0) * nLoads012[i][1](2)) / (sRef * cRef);
+        }
+        b->setUnsteadyMomentCoeff(imd, ifq, bcm1);
+
+        // compute load coefficients, only valid for small AoAs
+        sdpmVector3cd alpha(b->getAmplitude(imd), _pbl.getAoa(), b->getAmplitude(imd));
+        sdpmVectorXcd cfx(5);
+        cfx << cf[2].conjugate()(0), cf[1].conjugate()(0), cf[0](0), cf[1](0), cf[2](0);
+        sdpmVectorXcd cfz(5);
+        cfz << cf[2].conjugate()(2), cf[1].conjugate()(2), cf[0](2), cf[1](2), cf[2](2);
+        b->setUnsteadyDragCoeff(imd, ifq, (cf[1](0) * ca + cf[1](1) * sa + convolve(cfz, alpha)(4) * ca) / sRef);
+        b->setUnsteadySideforceCoeff(imd, ifq, (-cf[1](0) * sa + cf[1](1) * ca - convolve(cfz, alpha)(4) * sa) / sRef);
+        b->setUnsteadyLiftCoeff(imd, ifq, (-convolve(cfx, alpha)(4) + cf[1](2)) / sRef);
+
+        // compute total
+        _cd1[imd][ifq] += b->getUnsteadyDragCoeff(imd, ifq);
+        _cs1[imd][ifq] += b->getUnsteadySideforceCoeff(imd, ifq);
+        _cl1[imd][ifq] += b->getUnsteadyLiftCoeff(imd, ifq);
+        _cm1[imd][ifq] += b->getUnsteadyMomentCoeff(imd, ifq);
+    }
+}
+
+void Solver::write(std::ostream &out) const
+{
+    out << "sdpm::Solver with\n"
+        << _pbl;
+}