Version 0.1.0 (setup)

Version 0.1.0

Context

This MR setups the project and provides the basic code architecture, as well as test cases.
Two things remain to be done to have a minimal working release:

1. implement the computation of the aerodynamic influence coefficients
2. implement the import of the nonlinear solution and its interpolation on the Cartesian grid

Tests

Passed on ubuntu20.04 (python 3.8.2) and msys2 (python 3.8.1). Note: timers test might fail on msys

fpm/src/fSolver.cpp

+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "fSolver.h"
+#include "fBuilder.h"
+#include "fProblem.h"
+#include "fBody.h"
+
+#include "wMshData.h"
+#include "wNode.h"
+#include "wElement.h"
+#include "wMem.h"
+#include "wTag.h"
+#include "wResults.h"
+#include "wMshExport.h"
+#include "wCache.h"
+
+using namespace tbox;
+using namespace fpm;
+
+#define ANSI_COLOR_YELLOW "\x1b[1;33m"
+#define ANSI_COLOR_RESET "\x1b[0m"
+
+/**
+ * @brief Initialize the solver and perform sanity checks
+ * @authors Adrien Crovato
+ */
+Solver::Solver(std::shared_ptr<Builder> _aic) : aic(_aic), Cl(0), Cd(0), Cs(0), Cm(0)
+{
+    // Say hi
+    std::cout << std::endl;
+    std::cout << "*******************************************************************************" << std::endl;
+    std::cout << "**                                    \\_/                                    **" << std::endl;
+    std::cout << "**                           \\_______O(_)O_______/                           **" << std::endl;
+    std::cout << "**                                                                           **" << std::endl;
+    std::cout << "*******************************************************************************" << std::endl;
+    std::cout << "** Hi! My name is FPM v0.1-2007                                              **" << std::endl;
+    std::cout << "** Adrien Crovato & Romain Boman                                             **" << std::endl;
+    std::cout << "** ULiege 2020-2021                                                          **" << std::endl;
+    std::cout << "*******************************************************************************" << std::endl
+              << std::endl;
+
+    // Setup variables
+    phi.resize(aic->pbl->msh->nodes.size(), 0.);
+    vPhi.resize(aic->pbl->msh->nodes.size(), 0.);
+    U.resize(aic->pbl->msh->elems.size(), Eigen::Vector3d::Zero());
+    rho.resize(aic->pbl->msh->elems.size(), 0.);
+    M.resize(aic->pbl->msh->elems.size(), 0.);
+    Cp.resize(aic->pbl->msh->elems.size(), 0.);
+}
+
+/**
+ * @brief Solve the (full) potential equation
+ */
+void Solver::run()
+{
+    if (!aic->valid)
+        throw std::runtime_error("fpm::Solver:run() AIC matrices are not up-to-date!\n");
+
+    // @todo compute sigma from an interpolator and compute volume flow (interpolation should be performed only once BEFORE run)
+    std::cout << "Computing field sources... " << std::flush;
+    Eigen::VectorXd sigma(aic->nF);
+    if (aic->pbl->field)
+        sigma = Eigen::VectorXd::Zero(aic->nF);
+    Eigen::VectorXd uX = aic->Cx * sigma;
+    Eigen::VectorXd uY = aic->Cy * sigma;
+    Eigen::VectorXd uZ = aic->Cz * sigma;
+    std::cout << "done!" << std::endl;
+
+    std::cout << "Computing body sources... " << std::flush;
+    Eigen::VectorXd mu(aic->nP), tau(aic->nP);
+    int ig = 0;
+    for (auto body : aic->pbl->bodies)
+        for (int i = 0; i < body->tag->elems.size(); ++i, ++ig)
+            tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->normal());
+    std::cout << "done!" << std::endl;
+
+    std::cout << "Solving for body doublets... " << std::flush;
+    Eigen::PartialPivLU<Eigen::MatrixXd>(aic->A).solve(aic->B * tau);
+    std::cout << "done!" << std::endl;
+
+    std::cout << "Computing flow and loads on bodies... " << std::flush;
+    computeFlow(mu, tau, sigma);
+    computeLoad();
+    std::cout << "done!" << std::endl
+              << std::endl;
+}
+
+/**
+ * @brief Write the results
+ */
+void Solver::save(std::shared_ptr<MshExport> mshWriter, int n)
+{
+    // Write files
+    std::cout << "Saving files... " << std::endl;
+    // setup results
+    Results results;
+    results.scalars_at_nodes["phi"] = &phi;
+    results.scalars_at_nodes["varPhi"] = &vPhi;
+    results.vectors_at_elems["U"] = &U;
+    results.scalars_at_elems["rho"] = &rho;
+    results.scalars_at_elems["Mach"] = &M;
+    results.scalars_at_elems["Cp"] = &Cp;
+    // save (mesh and boundary surface)
+    if (n > 0)
+    {
+        mshWriter->save(aic->pbl->msh->name + "_" + std::to_string(n), results);
+        for (auto body : aic->pbl->bodies)
+            body->save(body->tag->name + "_" + std::to_string(n));
+    }
+    else
+    {
+        mshWriter->save(aic->pbl->msh->name, results);
+        for (auto body : aic->pbl->bodies)
+            body->save(body->tag->name);
+    }
+}
+
+/**
+ * @brief Compute flow solution on bodies
+ */
+void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, const Eigen::VectorXd &sigma)
+{
+    // @todo should we
+    // - compute velocity from phi at nodes or mu at nodes? if mu, then store phi at elements only
+
+    // Compute phi at elements and transfer at nodes
+    Eigen::VectorXd phie = aic->A * mu + aic->B * tau + aic->C * sigma;
+    auto pbl = aic->pbl;
+    // reset all values
+    for (auto n : pbl->msh->nodes)
+        vPhi[n->no - 1] = 0;
+    int ig = 0;
+    for (auto body : pbl->bodies)
+    {
+        // associate the potential value to its element
+        std::map<Element *, double> mapPhi;
+        for (auto e : body->tag->elems)
+        {
+            mapPhi[e] = phie[ig];
+            ig++;
+        }
+        // average at nodes
+        for (auto n : body->nodes)
+            for (auto e : body->neMap.at(n))
+                vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); // maybe use area average?
+    }
+
+    // Compute surface flow
+    for (auto body : pbl->bodies)
+    {
+        for (auto n : body->nodes)
+            phi[n->no - 1] = pbl->PhiI(n->pos) + vPhi[n->row]; // total potential
+
+        for (auto e : body->tag->elems)
+        {
+            U[e->no - 1] = pbl->U(*e, phi);     // velocity
+            rho[e->no - 1] = pbl->Rho(*e, phi); // density
+            M[e->no - 1] = pbl->Mach(*e, phi);  // Mach number
+            Cp[e->no - 1] = pbl->Cp(*e, phi);   // pressure coefficient
+        }
+    }
+}
+
+/**
+ * @brief Compute aerodynamic load coefficients on bodies
+ */
+void Solver::computeLoad()
+{
+    Cl = 0;
+    Cd = 0;
+    Cs = 0;
+    Cm = 0;
+    auto pbl = aic->pbl;
+    for (auto body : pbl->bodies)
+    {
+        // Compute normalized load (i.e. load / pressure)
+        for (size_t i = 0; i < body->nodes.size(); ++i)
+        {
+            body->cLoadX[i] = 0;
+            body->cLoadY[i] = 0;
+            body->cLoadZ[i] = 0;
+            for (auto e : body->neMap.at(body->nodes[i]))
+            {
+                Eigen::Vector3d cLoad = Cp[e->no - 1] * e->getVMem().getVol() * e->normal() / e->nodes.size();
+                body->cLoadX[i] += cLoad(0);
+                body->cLoadY[i] += cLoad(1);
+                body->cLoadZ[i] += cLoad(2);
+            }
+        }
+
+        // Compute aerodynamic load coefficients (i.e. Load / pressure / surface)
+        // compute coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
+        body->Cm = 0;
+        Eigen::Vector3d Cf(0, 0, 0);
+        for (auto e : body->tag->elems)
+        {
+            Eigen::Vector3d l = e->getVMem().getCg() - pbl->xR;                       // lever arm
+            Eigen::Vector3d cf = Cp[e->no - 1] * e->getVMem().getVol() * e->normal(); // force coefficient
+            Cf -= cf;
+            body->Cm -= (l(2) * cf(0) - l(0) * cf(2)) / (pbl->sR * pbl->cR); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+        }
+        // rotate to flow direction
+        body->Cd = (Cf(0) * cos(pbl->alpha) * cos(pbl->beta) + Cf(1) * sin(pbl->beta) + Cf(2) * sin(pbl->alpha) * cos(pbl->beta)) / pbl->sR;
+        body->Cs = (-Cf(0) * cos(pbl->alpha) * sin(pbl->beta) + Cf(1) * cos(pbl->beta) - Cf(2) * sin(pbl->alpha) * sin(pbl->beta)) / pbl->sR;
+        body->Cl = (-Cf(0) * sin(pbl->alpha) + Cf(2) * cos(pbl->alpha)) / pbl->sR;
+        // compute total
+        Cl += body->Cl;
+        Cd += body->Cd;
+        Cs += body->Cs;
+        Cm += body->Cm;
+    }
+}
+
+void Solver::write(std::ostream &out) const
+{
+    out << "fpm::Solver\nwith\t" << *aic;
+}