From 4981213b5fa06cd4cc0bbdea46ea3410a0b8fd3e Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Sat, 17 Feb 2024 12:30:58 +0100
Subject: [PATCH] (feat) First part for coupled adjoint VII
---
blast/_src/blastw.h | 3 +-
blast/_src/blastw.i | 5 +-
blast/interfaces/dart/DartInterface.py | 10 +-
blast/src/DCoupledAdjoint.cpp | 277 +++++++++++++++++++++++++
blast/src/DCoupledAdjoint.h | 94 +++++++++
blast/src/blast.h | 3 +
blast/tests/dart/adjointVII.py | 218 +++++++++++++++++++
7 files changed, 606 insertions(+), 4 deletions(-)
create mode 100644 blast/src/DCoupledAdjoint.cpp
create mode 100644 blast/src/DCoupledAdjoint.h
create mode 100644 blast/tests/dart/adjointVII.py
diff --git a/blast/_src/blastw.h b/blast/_src/blastw.h
index ec45463..b1c27f1 100644
--- a/blast/_src/blastw.h
+++ b/blast/_src/blastw.h
@@ -14,4 +14,5 @@
* limitations under the License.
*/
-#include "DDriver.h"
\ No newline at end of file
+#include "DDriver.h"
+#include "DCoupledAdjoint.h"
\ No newline at end of file
diff --git a/blast/_src/blastw.i b/blast/_src/blastw.i
index ab103f1..ac1b2e6 100644
--- a/blast/_src/blastw.i
+++ b/blast/_src/blastw.i
@@ -54,4 +54,7 @@ threads="1"
%immutable blast::Driver::Cdp_sec;
%immutable blast::Driver::CFL0;
%immutable blast::Driver::verbose;
-%include "DDriver.h"
\ No newline at end of file
+%include "DDriver.h"
+
+%shared_ptr(blast::CoupledAdjoint);
+%include "DCoupledAdjoint.h"
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index da0809b..a52b126 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -28,10 +28,16 @@ class DartInterface(SolversInterface):
def __init__(self, dartCfg, vSolver, _cfg):
try:
from modules.dartflo.dart.api.core import initDart
- self.argOut = initDart(dartCfg, viscous=True)
+ if 'task' not in dartCfg:
+ dartCfg['task'] = 'analysis'
+ self.argOut = initDart(dartCfg, task=dartCfg['task'], viscous=True)
self.solver = self.argOut.get('sol') # Solver
self.blw = [self.argOut.get('blwb'), self.argOut.get('blww')]
- print('Dart successfully loaded.')
+ if dartCfg['task'] == 'optimization':
+ self.adjointSolver = self.argOut.get('adj')
+ print('Dart successfully loaded in optimization mode.')
+ else:
+ print('Dart successfully loaded in analysis mode.')
except:
print(ccolors.ANSI_RED, 'Could not load DART. Make sure it is installed.', ccolors.ANSI_RESET)
raise RuntimeError('Import error')
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
new file mode 100644
index 0000000..34db14b
--- /dev/null
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -0,0 +1,277 @@
+/*
+ * 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 "DCoupledAdjoint.h"
+#include "../../modules/dartflo/dart/src/wAdjoint.h"
+#include "../../modules/dartflo/dart/src/wNewton.h"
+#include "../../modules/dartflo/dart/src/wSolver.h"
+#include "../../modules/dartflo/dart/src/wProblem.h"
+#include "../../modules/dartflo/dart/src/wBlowing.h"
+#include "../../modules/dartflo/dart/src/wFluid.h"
+#include "wMshData.h"
+#include "wNode.h"
+
+#include "wLinearSolver.h"
+#include "wGmres.h"
+
+#include <Eigen/Sparse>
+#include <deque>
+#include <algorithm>
+
+using namespace blast;
+
+CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::shared_ptr<blast::Driver> vSolver) : iadjoint(_iadjoint), vsol(vSolver)
+{
+ isol = iadjoint->sol;
+ // Linear solvers (if GMRES, use the same, but with a thighter tolerance)
+ if (isol->linsol->type() == tbox::LSolType::GMRES_ILUT)
+ {
+ std::shared_ptr<tbox::Gmres> gmres = std::make_shared<tbox::Gmres>(*dynamic_cast<tbox::Gmres *>(isol->linsol.get()));
+ gmres->setTolerance(1e-8);
+ alinsol = gmres;
+ }
+ else
+ alinsol = isol->linsol;
+
+ sizeBlowing = iadjoint->getSizeBlowing();
+
+ dCldAoA = 0.0;
+ dCddAoA = 0.0;
+}
+
+void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint)
+{
+ // Determine the sizes of the smaller matrices
+ int rows = matrixAdjoint.rows();
+ int cols = matrixAdjoint.cols();
+
+ // Create a list of triplets for the larger matrix
+ std::vector<Eigen::Triplet<double>> triplets;
+
+ // Iterate over the rows and columns of the vector
+ int rowOffset = 0;
+ for (const auto& row : matrices) {
+ int colOffset = 0;
+ for (const auto& matrix : row) {
+ // Add the triplets of the matrix to the list of triplets for the larger matrix
+ for (int k = 0; k < matrix->outerSize(); ++k) {
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(*matrix, k); it; ++it) {
+ triplets.push_back(Eigen::Triplet<double>(it.row() + rowOffset, it.col() + colOffset, it.value()));
+ }
+ }
+ colOffset += matrix->cols();
+ }
+ rowOffset += row[0]->rows();
+ }
+
+ // Create a new matrix from the deque of triplets
+ matrixAdjoint.setFromTriplets(triplets.begin(), triplets.end());
+}
+
+void CoupledAdjoint::run()
+{
+ Eigen::SparseMatrix<double, Eigen::RowMajor> A = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> B = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> C = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> D = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> E = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> F = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> G = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> H = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> I = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
+
+ std::deque<Eigen::Triplet<double>> T;
+ T.push_back(Eigen::Triplet<double>(0, 0, 1.0));
+ T.push_back(Eigen::Triplet<double>(1, 1, 2.0));
+ A.setFromTriplets(T.begin(), T.end());
+ //std::cout << "A\n" <<A << std::endl;
+ T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 3.0));
+ T.push_back(Eigen::Triplet<double>(0, 1, 2.0));
+ B.setFromTriplets(T.begin(), T.end());
+ //std::cout << "B\n" <<B << std::endl;
+ /*T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 0.0));
+ T.push_back(Eigen::Triplet<double>(1, 1, 0.0));
+ C.setFromTriplets(T.begin(), T.end());*/
+ std::cout << "C\n" << C << std::endl;
+ T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
+ T.push_back(Eigen::Triplet<double>(2, 1, 1.0));
+ D.setFromTriplets(T.begin(), T.end());
+ //std::cout << "D\n" << D << std::endl;
+
+ T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
+ T.push_back(Eigen::Triplet<double>(1, 1, 1.0));
+ E.setFromTriplets(T.begin(), T.end());
+ //std::cout << "E\n" << E << std::endl;
+
+ T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
+ T.push_back(Eigen::Triplet<double>(1, 1, 1.0));
+ F.setFromTriplets(T.begin(), T.end());
+ //std::cout << "F\n" << F << std::endl;
+
+ T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
+ T.push_back(Eigen::Triplet<double>(1, 1, 4.0));
+ G.setFromTriplets(T.begin(), T.end());
+ //std::cout << "G\n" << G << std::endl;
+
+ T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
+ T.push_back(Eigen::Triplet<double>(1, 1, 1.0));
+ H.setFromTriplets(T.begin(), T.end());
+ //std::cout << "H\n" << H << std::endl;
+
+ T.clear();
+ T.push_back(Eigen::Triplet<double>(0, 0, 5.0));
+ T.push_back(Eigen::Triplet<double>(1, 2, 1.0));
+ I.setFromTriplets(T.begin(), T.end());
+ //std::cout << "I\n" << I << std::endl;
+
+ std::cout << "gradientswrtInviscidFlow\n";
+ dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->nodes.size(), isol->pbl->msh->nodes.size());
+ dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->nodes.size(), isol->pbl->msh->nodes.size());
+ gradientswrtInviscidFlow();
+ std::cout << " coucou " << std::endl;
+ std::cout << "dRM_phi\n" << dRM_phi.norm() << std::endl;
+ std::cout << "dRrho_phi\n" << dRrho_phi.norm() << std::endl;
+
+ // Store pointers to the matrices in a 2D vector
+ std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> matrices = {
+ {&A, &B, &C, &D, &E},
+ {&D, &E, &F, &I, &G},
+ {&G, &H, &I, &A, &B},
+ {&B, &E, &F, &I, &G},
+ {&B, &A, &I, &C, &B}
+ };
+
+ // Create a larger sparse matrix and set it from the list of triplets
+ int rows = 0; int cols = 0;
+ for (const auto& row : matrices) {
+ rows += row[0]->rows(); // All matrices in a row have the same number of rows
+ cols = std::max(cols, static_cast<int>(row.size() * row[0]->cols())); // All matrices in a column have the same number of columns
+ }
+ Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols);
+
+ buildAdjointMatrix(matrices, A_adjoint);
+
+ std::cout << "largeMatrix\n" << A_adjoint << std::endl;
+
+ // Solve coupled system
+ // alinsol->compute(AdjointMatrix.transpose(), Eigen::Map<Eigen::VectorXd>(dU_.data(), dU_.size()), lambdas);
+}
+
+void CoupledAdjoint::gradientswrtInviscidFlow()
+{
+ //**** dRphi_phi ****//
+ dRphi_phi = iadjoint->getRu_U();
+
+ //**** dRM_phi, dRrho_phi, dRv_phi ****//
+ auto pbl = isol->pbl;
+ //dRM_dPhi = iadjoint->getRM_U();
+ // Finite difference
+ double deltaPhi = 1e-4; // perturbation
+ std::vector<double> Phi_save = isol->phi; // Differential of the independant variable (phi)
+ std::vector<double> dM = std::vector<double>(pbl->bBCs[0]->nodes.size(), 0.0); // Differential of the Mach number
+ std::vector<double> drho = std::vector<double>(pbl->bBCs[0]->nodes.size(), 0.0); // Differential of the density
+
+ std::vector<Eigen::Triplet<double>> tripletsdM;
+ std::vector<Eigen::Triplet<double>> tripletsdrho;
+
+ for (auto n : pbl->msh->nodes){
+ Phi_save[n->row] = isol->phi[n->row] + deltaPhi;
+
+ std::fill(dM.begin(), dM.end(), 0.0);
+ std::fill(drho.begin(), drho.end(), 0.0);
+
+ // Recompute Mach number on surface nodes.
+ size_t jj = 0;
+ for (auto srfNode : pbl->bBCs[0]->nodes){
+ // Average of the Mach on the elements adjacent to the considered node.
+ for (auto e : pbl->fluid->neMap[srfNode]){
+ dM[jj] += pbl->fluid->mach->eval(*e, Phi_save, 0);
+ drho[jj] += pbl->fluid->rho->eval(*e, Phi_save, 0);
+ }
+ dM[jj] /= pbl->fluid->neMap[srfNode].size();
+ drho[jj] /= pbl->fluid->neMap[srfNode].size();
+ // Form triplets
+ tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row, (dM[jj] - isol->M[srfNode->row])/deltaPhi));
+ tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row, (drho[jj] - isol->rho[srfNode->row])/deltaPhi));
+ jj++;
+ }
+ // Reset phi
+ Phi_save[n->row] = isol->phi[n->row];
+ }
+
+ // Fill matrices with gradients
+ dRM_phi.setFromTriplets(tripletsdM.begin(), tripletsdM.end());
+ dRrho_phi.setFromTriplets(tripletsdrho.begin(), tripletsdrho.end());
+ // Remove zeros
+ dRM_phi.prune(0.);
+ dRrho_phi.prune(0.);
+
+ // // Objective functionss
+ // dCl_dPhi =
+ // dCd_dPhi =
+}
+
+void CoupledAdjoint::gradientswrtMachNumber(){
+ // dRphi_dM = 0
+ // dRM_dM =
+ // dRrho_dM = 0
+ // dRv_dM = 0
+ // dRblw_dM =
+
+ // // Objective functions
+ // dCl_dM =
+ // dCd_dM =
+}
+
+void CoupledAdjoint::gradientswrtDensity(){
+ // dRphi_dRho = 0
+ // dRM_dRho = 0
+ // dRrho_dRho =
+ // dRv_dRho = 0
+ // dRblw_dRho =
+
+ // dCl_dRho =
+ // dCd_dRho =
+}
+
+void CoupledAdjoint::gradientswrtVelocity(){
+ // dRphi_dV = 0
+ // dRM_dV = 0
+ // dRrho_dV = 0
+ // dRv_dV =
+ // dRblw_dV =
+
+ // dCl_dV =
+ // dCd_dV =
+}
+
+void CoupledAdjoint::gradientswrtBlowingVelocity(){
+ // dRphi_dBlw =
+ // dRM_dBlw = 0
+ // dRrho_dBlw = 0
+ // dRv_dBlw = 0
+ // dRblw_dBlw =
+
+ // dCl_dBlw =
+ // dCd_dBlw =
+}
\ No newline at end of file
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
new file mode 100644
index 0000000..73aa0b5
--- /dev/null
+++ b/blast/src/DCoupledAdjoint.h
@@ -0,0 +1,94 @@
+/*
+ * 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.
+ */
+
+#ifndef DCOUPLEDADJOINT_H
+#define DCOUPLEDADJOINT_H
+
+#include "blast.h"
+#include "wObject.h"
+#include "../../modules/dartflo/dart/src/wAdjoint.h"
+#include <Eigen/Sparse>
+
+#include <iostream>
+
+namespace blast
+{
+
+/**
+ * @brief Adjoint of the coupled system.
+ * @author Paul Dechamps
+ */
+class BLAST_API CoupledAdjoint : public fwk::wSharedObject
+{
+
+public:
+ std::shared_ptr<dart::Adjoint> iadjoint; ///< Adjoint solver
+ std::shared_ptr<dart::Newton> isol; ///< Inviscid Newton solver
+ std::shared_ptr<blast::Driver> vsol; ///< Viscous solver
+ std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero)
+ double dCldAoA;
+ double dCddAoA;
+
+private:
+ size_t sizeBlowing;
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_phi; ///< Inviscid flow Jacobian
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_M; ///< Derivative of the inviscid flow residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_v; ///< Derivative of the inviscid flow residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_rho; ///< Derivative of the inviscid flow residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_blw; ///< Derivative of the inviscid flow residual with respect to the blowing velocity
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi; ///< Derivative of the Mach number residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_M; ///< Derivative of the Mach number residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_v; ///< Derivative of the Mach number residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_rho; ///< Derivative of the Mach number residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_blw; ///< Derivative of the Mach number residual with respect to the blowing velocity
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_phi; ///< Derivative of the density residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_M; ///< Derivative of the density residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_v; ///< Derivative of the density residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_rho; ///< Derivative of the density residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_blw; ///< Derivative of the density residual with respect to the blowing velocity
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_phi; ///< Derivative of the velocity residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_M; ///< Derivative of the velocity residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_v; ///< Derivative of the velocity residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_rho; ///< Derivative of the velocity residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_blw; ///< Derivative of the velocity residual with respect to the blowing velocity
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_phi; ///< Derivative of the blowing velocity residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_M; ///< Derivative of the blowing velocity residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_v; ///< Derivative of the blowing velocity residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_rho; ///< Derivative of the blowing velocity residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_blw; ///< Derivative of the blowing velocity residual with respect to the blowing velocity
+
+public:
+ CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
+ virtual ~CoupledAdjoint () {std::cout << "~CoupledAdjoint()\n";};
+
+ void run();
+
+private:
+ void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
+
+ void gradientswrtInviscidFlow();
+ void gradientswrtMachNumber();
+ void gradientswrtDensity();
+ void gradientswrtVelocity();
+ void gradientswrtBlowingVelocity();
+};
+} // namespace blast
+#endif // DDARTADJOINT
diff --git a/blast/src/blast.h b/blast/src/blast.h
index 71cdab8..6b166f5 100644
--- a/blast/src/blast.h
+++ b/blast/src/blast.h
@@ -48,6 +48,9 @@ class Edge;
// Other
class Closures;
+
+// Adjoint
+class CoupledAdjoint;
} // namespace blast
#endif // BLAST_H
diff --git a/blast/tests/dart/adjointVII.py b/blast/tests/dart/adjointVII.py
new file mode 100644
index 0000000..6b2f8bb
--- /dev/null
+++ b/blast/tests/dart/adjointVII.py
@@ -0,0 +1,218 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2022 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.
+
+
+# @author Paul Dechamps
+# @date 2022
+# Test the blast implementation. The test case is a compressible attached transitional flow.
+# Tested functionalities;
+# - Time integration.
+# - Two time-marching techniques (agressive and safe).
+# - Transition routines.
+# - Forced transition.
+# - Compressible flow routines.
+# - Laminar and turbulent flow.
+#
+# Untested functionalities.
+# - Separation routines.
+# - Multiple failure case at one iteration.
+# - Response to unconverged inviscid solver.
+
+# Imports.
+
+import blast.utils as viscUtils
+import blast
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+ import os.path
+ # Parameters
+ return {
+ # Options
+ 'scenario' : 'aerodynamic',
+ 'task' : 'analysis',
+ 'Threads' : nthrds, # number of threads
+ 'Verb' : verb, # verbosity
+ # Model (geometry or mesh)
+ 'File' : os.path.dirname(os.path.abspath(__file__)) + '/../../models/dart/n0012.geo', # Input file containing the model
+ 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 2, 'msTe' : 0.05, 'msLe' : 0.05}, # parameters for input file model
+ 'Dim' : 2, # problem dimension
+ 'Format' : 'gmsh', # save format (vtk or gmsh)
+ # Markers
+ 'Fluid' : 'field', # name of physical group containing the fluid
+ 'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+ 'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+ 'Wake' : 'wake', # LIST of names of physical group containing the wake
+ 'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+ 'Te' : 'te', # LIST of names of physical group containing the trailing edge
+ 'dbc' : True,
+ 'Upstream' : 'upstream',
+ # Freestream
+ 'M_inf' : 0.2, # freestream Mach number
+ 'AoA' : 2., # freestream angle of attack
+ # Geometry
+ 'S_ref' : 1., # reference surface length
+ 'c_ref' : 1., # reference chord length
+ 'x_ref' : .25, # reference point for moment computation (x)
+ 'y_ref' : 0.0, # reference point for moment computation (y)
+ 'z_ref' : 0.0, # reference point for moment computation (z)
+ # Numerical
+ 'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+ 'G_fill' : 2, # fill-in factor for GMRES preconditioner
+ 'G_tol' : 1e-5, # tolerance for GMRES
+ 'G_restart' : 50, # restart for GMRES
+ 'Rel_tol' : 1e-6, # relative tolerance on solver residual
+ 'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+ 'Max_it' : 20, # solver maximum number of iterations
+
+ 'task': 'optimization'
+ }
+
+def cfgBlast(verb):
+ return {
+ 'Re' : 1e7, # Freestream Reynolds number
+ 'Minf' : 0.2, # Freestream Mach number (used for the computation of the time step only)
+ 'CFL0' : 1, # Inital CFL number of the calculation
+ 'Verb': verb, # Verbosity level of the solver
+ 'couplIter': 25, # Maximum number of iterations
+ 'couplTol' : 1e-4, # Tolerance of the VII methodology
+ 'iterPrint': 5, # int, number of iterations between outputs
+ 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
+ 'sections' : [0],
+ 'xtrF' : [None, 0.4],# Forced transition location
+ 'interpolator' : 'Matching',
+ 'nDim' : 2
+ }
+
+def main():
+ # Timer.
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ tms['pre'].start()
+ coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ tms['pre'].stop()
+
+ adjointSolver = blast.CoupledAdjoint(iSolverAPI.adjointSolver, vSolver)
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ aeroCoeffs = coupler.run()
+ #iSolverAPI.run()
+
+ adjointSolver.run()
+ print("sucess", adjointSolver)
+
+ quit()
+
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ tms['pre'].start()
+ coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ tms['pre'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ #aeroCoeffs = coupler.run()
+ iSolverAPI.run()
+ quit()
+ tms['solver'].stop()
+
+ # Display results.
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' Re M alpha Cl Cd Cdp Cdf Cm')
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, iSolverAPI.getMinf(), iSolverAPI.getAoA()*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
+
+ # Write results to file.
+ vSolution = viscUtils.getSolution(vSolver)
+ vSolution['Cdt_int'] = vSolution['Cdf'] + iSolverAPI.getCd()
+ tms['total'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+ print('SOLVERS statistics')
+ print(coupler.tms)
+
+ # Test solution
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('Cl', iSolverAPI.getCl(), 0.230, 5e-2)) # XFOIL 0.2325
+ tests.add(CTest('Cd wake', vSolution['Cdt_w'], 0.0056, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cd integral', vSolution['Cdt_int'], 0.0059, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cdf', vSolution['Cdf'], 0.00465, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
+ tests.add(CTest('xtr_top', vSolution['xtrT'], 0.20, 0.03, forceabs=True)) # XFOIL 0.1877
+ tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
+ tests.add(CTest('Iterations', len(aeroCoeffs), 17, 0, forceabs=True))
+ tests.run()
+
+ # Show results
+ if not args.nogui:
+ iCp = viscUtils.read('Cp_inviscid.dat')
+ vCp = viscUtils.read('Cp_viscous.dat')
+ xfoilCp = viscUtils.read('../../blast/models/references/cpXfoil_n0012.dat', delim=None, skip = 1)
+ xfoilCpInv = viscUtils.read('../../blast/models/references/cpXfoilInv_n0012.dat', delim=None, skip = 1)
+ plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
+ np.column_stack((xfoilCp[:,0], xfoilCp[:,1])),
+ np.column_stack((iCp[:,0], iCp[:,3])),
+ np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
+ 'labels': ['Blast (VII)', 'XFOIL VII', 'DART (inviscid)', 'XFOIL (inviscid)'],
+ 'lw': [3, 3, 2, 2],
+ 'color': ['firebrick', 'darkblue', 'firebrick', 'darkblue'],
+ 'ls': ['-', '-', '--', '--'],
+ 'reverse': True,
+ 'xlim':[0, 1],
+ 'yreverse': True,
+ 'legend': True,
+ 'xlabel': '$x/c$',
+ 'ylabel': '$c_p$'
+ }
+ viscUtils.plot(plotcp)
+
+ xfoilBl = viscUtils.read('../../blast/models/references/blXfoil_n0012.dat', delim=None, skip = 1)
+ plotcf = {'curves': [np.column_stack((vSolution['x'], vSolution['Cf'])),
+ np.column_stack((xfoilBl[:,1], xfoilBl[:,6]))],
+ 'labels': ['Blast', 'XFOIL'],
+ 'lw': [3, 3],
+ 'color': ['firebrick', 'darkblue'],
+ 'ls': ['-', '-'],
+ 'reverse': True,
+ 'xlim':[0, 1],
+ 'legend': True,
+ 'xlabel': '$x/c$',
+ 'ylabel': '$c_f$'
+ }
+ viscUtils.plot(plotcf)
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
--
GitLab