Skip to content
Snippets Groups Projects
Verified Commit 4981213b authored by Paul Dechamps's avatar Paul Dechamps :speech_balloon:
Browse files

(feat) First part for coupled adjoint VII

parent 91b08b19
No related branches found
No related tags found
No related merge requests found
Pipeline #20721 failed
......@@ -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
......@@ -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"
......@@ -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')
......
/*
* 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
/*
* 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
......@@ -48,6 +48,9 @@ class Edge;
// Other
class Closures;
// Adjoint
class CoupledAdjoint;
} // namespace blast
#endif // BLAST_H
#!/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()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment