diff --git a/CMake/FindCODI.cmake b/CMake/FindCODI.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..4e3551cf73d49ee5e7929acf7eab4b0d0504c760
--- /dev/null
+++ b/CMake/FindCODI.cmake
@@ -0,0 +1,75 @@
+# 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.
+
+# Try to find CoDiPack headers
+# ----------------------------------------------------------------------------
+# output:
+# CODI_FOUND : whether the packae was found
+# CODI_INCLUDE_DIRS : where the headers are [cached]
+# CODI_VERSION : version number
+# ----------------------------------------------------------------------------
+# autodetection:
+# use "CMAKE_PREFIX_PATH=/path/to/package/"
+# or "CMAKE_INCLUDE_PATH=/path/to/package/include"
+# or "INCLUDE=/path/to/package/include"
+# ----------------------------------------------------------------------------
+
+# Check version
+macro(_check_version)
+ file(READ "${CODI_INCLUDE_DIRS}/codi.hpp" _version_header)
+
+ string(REGEX MATCH "define[ \t]+CODI_MAJOR_VERSION[ \t]+([0-9]+)" major_version_match "${_version_header}")
+ set(CODI_MAJOR_VERSION "${CMAKE_MATCH_1}")
+ string(REGEX MATCH "define[ \t]+CODI_MINOR_VERSION[ \t]+([0-9]+)" minor_version_match "${_version_header}")
+ set(CODI_MINOR_VERSION "${CMAKE_MATCH_1}")
+ string(REGEX MATCH "define[ \t]+CODI_BUILD_VERSION[ \t]+([0-9]+)" build_version_match "${_version_header}")
+ set(CODI_BUILD_VERSION "${CMAKE_MATCH_1}")
+
+ set(CODI_VERSION ${CODI_MAJOR_VERSION}.${CODI_MINOR_VERSION}.${CODI_BUILD_VERSION})
+
+ if(${CODI_VERSION} VERSION_LESS ${CODI_FIND_VERSION})
+ message(FATAL_ERROR "CoDiPack version ${CODI_VERSION} found in ${CODI_INCLUDE_DIRS}, but at least version ${CODI_FIND_VERSION} is required!")
+ endif()
+endmacro()
+
+# Set a dummy version if not provided
+if(NOT CODI_FIND_VERSION)
+ if(NOT CODI_FIND_VERSION_MAJOR)
+ set(CODI_FIND_VERSION_MAJOR 2)
+ endif()
+ if(NOT CODI_FIND_VERSION_MINOR)
+ set(CODI_FIND_VERSION_MINOR 0)
+ endif()
+ if(NOT CODI_FIND_VERSION_BUILD)
+ set(CODI_FIND_VERSION_BUILD 2)
+ endif()
+
+ set(CODI_FIND_VERSION "${CODI_FIND_VERSION_MAJOR}.${CODI_FIND_VERSION_MINOR}.${CODI_FIND_VERSION_BUILD}")
+endif()
+
+# Find the header and check the version
+find_path(CODI_INCLUDE_DIRS "codi.hpp")
+if (CODI_INCLUDE_DIRS)
+ _check_version()
+else()
+ message(FATAL_ERROR "CoDiPack headers not found! Define the path in the INCLUDE environment variable.")
+endif()
+
+# handle the QUIETLY and REQUIRED arguments and set EIGEN_FOUND to TRUE
+# if all listed variables are TRUE
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(CODI
+ FOUND_VAR CODI_FOUND
+ REQUIRED_VARS CODI_INCLUDE_DIRS
+ VERSION_VAR CODI_VERSION)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 03edcc4736eda7a6192931ff055e8bfb4f592861..b916eb53c6c6efddc4d4248a95ba804b9ff3c1b9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -54,6 +54,8 @@ ELSEIF(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
#SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Weverything") # add verbosity
ENDIF()
+OPTION(USE_CODI "Compile with CoDiPack" OFF)
+
# -- DEPENDENCIES
# Python
# use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
@@ -61,7 +63,14 @@ FIND_PACKAGE(Python3 COMPONENTS Interpreter Development)
SET(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
SET(PYTHON_LIBRARIES ${Python3_LIBRARIES})
SET(PYTHON_INCLUDE_PATH ${Python3_INCLUDE_DIRS})
-SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})
+SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})
+
+# CoDiPack
+IF(USE_CODI)
+ FIND_PACKAGE(CODI 2.1.0 REQUIRED)
+ MESSAGE(STATUS "CODI_INCLUDE_DIRS: ${CODI_INCLUDE_DIRS}")
+ENDIF()
+
# SWIG
FIND_PACKAGE(SWIG REQUIRED)
@@ -74,6 +83,10 @@ ENDIF()
# -- DEFINE (for SWIG to detect definitions)
INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}) # to find "amfe_def.h"
+# -- DEFINE (for SWIG to detect definitions)
+CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/blaster_def.h.in ${CMAKE_BINARY_DIR}/blaster_def.h)
+INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}) # pass to the most top-level build dir to find "sdpm_def.h"
+
# -- CTest
ENABLE_TESTING()
@@ -109,3 +122,4 @@ MESSAGE(STATUS "* BUILD TYPE: ${CMAKE_BUILD_TYPE}")
MESSAGE(STATUS "* VTK SUPPORT: ${USE_VTK}")
MESSAGE(STATUS "* MKL SUPPORT: ${USE_MKL}")
MESSAGE(STATUS "* MUMPS SUPPORT: ${USE_MUMPS}")
+MESSAGE(STATUS "* CoDiPack SUPPORT: ${USE_CODI}")
diff --git a/blast/_src/blastw.i b/blast/_src/blastw.i
index 7938cd7f1ffda01b03e9a21d929dceb40670685c..58432beb1fc04e8bc5faf87b4b5cbe1379be8bb2 100644
--- a/blast/_src/blastw.i
+++ b/blast/_src/blastw.i
@@ -37,6 +37,39 @@ threads="1"
%}
+%include "std_vector.i"
+
+%include "blaster_def.h" // needed to access USE_CODI in SWIG interface
+#ifndef USE_CODI
+typedef double bldouble; // needed so that SWIG does not convert bldouble to double
+#endif
+%typemap(in) bldouble* (bldouble temp) {
+ if (PyFloat_Check($input)) {
+ temp = bldouble(PyFloat_AsDouble($input));
+ $1 = &temp;
+ } else if (PyObject_TypeCheck($input, &PyCapsule_Type)) {
+ $1 = (bldouble *) PyCapsule_GetPointer($input, NULL);
+ } else if (SWIG_IsOK(SWIG_ConvertPtr($input, (void **)&$1, SWIGTYPE_p_bldouble, 0))) {
+ // Successfully converted to bldouble*
+ } else {
+ PyErr_SetString(PyExc_TypeError, "Expected a float or bldouble");
+ return NULL;
+ }
+}
+%typemap(out) bldouble {
+#ifndef USE_CODI
+ $result = PyFloat_FromDouble($1);
+#else
+ $result = PyFloat_FromDouble(($1).getValue());
+#endif
+}
+#ifndef USE_CODI
+%rename(std_vector_double) std::vector<bldouble>;
+%rename(std_vector_vector_double) std::vector<std::vector<bldouble>>;
+#else
+%template(std_vector_bldouble) std::vector<bldouble>;
+%template(std_vector_vector_bldouble) std::vector<std::vector<bldouble>>;
+#endif
%include "fwkw.swg"
@@ -78,6 +111,74 @@ threads="1"
%template(std_vector_blSection) std::vector<std::shared_ptr<blast::Section>>;
%template(std_vector_blBoundaryLayer) std::vector<std::shared_ptr<blast::BoundaryLayer>>;
+#ifdef USE_CODI
+class bldouble {};
+%extend bldouble {
+ double __float__()
+ {
+ return (*self).getValue();
+ }
+ bldouble __add__(bldouble const &other)
+ {
+ return *self + other;
+ }
+ bldouble __sub__(bldouble const &other)
+ {
+ return *self - other;
+ }
+ bldouble __mul__(bldouble const &other)
+ {
+ return *self * other;
+ }
+ bldouble __div__(bldouble const &other)
+ {
+ return *self / other;
+ }
+ // use repr instead of str so that iterables also print nicely
+ std::string __repr__()
+ {
+ std::ostringstream out; out << *self;
+ return out.str();
+ }
+}
+#endif
+
+class blVector3d {};
+%extend blVector3d {
+ blVector3d(double x, double y, double z)
+ {
+ blVector3d *v = new blVector3d();
+ *v << x,y,z;
+ return v;
+ }
+ double __getitem__(int i)
+ {
+ if (i > 2)
+ throw std::runtime_error("blVector3d: index out of range (size 3)!\n");
+#ifndef USE_CODI
+ return (*self)(i);
+#else
+ return (*self)(i).getValue();
+#endif
+ }
+ // use repr instead of str so that iterables also print nicely
+ std::string __repr__()
+ {
+ std::ostringstream out; out << *self;
+ return out.str();
+ }
+}
+
+%extend std::vector<bldouble> {
+ void __setitem__(size_t i, double val) {
+#ifndef USE_CODI
+ (*$self)[i] = val;
+#else
+ (*$self)[i] = static_cast<bldouble>(val); // Explicit conversion
+#endif
+ }
+}
+
%extend blast::Section {
%pythoncode {
def __iter__(self):
@@ -92,4 +193,4 @@ threads="1"
else:
raise StopIteration
}
-}
\ No newline at end of file
+}
diff --git a/blast/blCoupler.py b/blast/blCoupler.py
index e03334f74613283abe2e922fe5a64cfcc57c9dfd..8810afbcda88b786f2731b785c84c7a0e8f17cf0 100644
--- a/blast/blCoupler.py
+++ b/blast/blCoupler.py
@@ -113,12 +113,12 @@ class Coupler:
self.isol.getBlowingBoundary(couplIter+1)
aeroCoeffs['Cl'].append(self.isol.getCl())
- aeroCoeffs['Cd'].append(self.isol.getCd() + self.vsol.Cdf)
- aeroCoeffs['Cdwake'].append(self.vsol.Cdt)
+ aeroCoeffs['Cd'].append(self.isol.getCd() + self.vsol.getCdf())
+ aeroCoeffs['Cdwake'].append(self.vsol.getCdt())
# Check convergence status.
- cd = self.vsol.Cdf + self.isol.getCd()
- #cd = self.vsol.Cdt if self.vsol.Cdt != 0 else self.vsol.Cdf + self.isol.getCd()
+ cd = self.vsol.getCdf() + self.isol.getCd()
+ #cd = self.vsol.getCdt() if self.vsol.getCdt() != 0 else self.vsol.getCdf() + self.isol.getCd()
error = abs((cd - cdPrev) / cd) if cd != 0 else 1
# Save restart file
@@ -127,7 +127,7 @@ class Coupler:
self.isol.save_restart(cd, sfx=_sfx)
if error <= self.tol:
- print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>5.0f} {:>5.0f} | {:>6.3f}\n'.format(couplIter, self.isol.getCl(), self.isol.getTotalDrag(), self.vsol.Cdt, self.vsol.bodies[0].getAvgxtr(0), self.vsol.bodies[0].getAvgxtr(1), iEc, vEc, np.log10(error)), ccolors.ANSI_RESET)
+ print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>5.0f} {:>5.0f} | {:>6.3f}\n'.format(couplIter, self.isol.getCl(), self.isol.getTotalDrag(), self.vsol.getCdt(), self.vsol.bodies[0].getAvgxtr(0), self.vsol.bodies[0].getAvgxtr(1), iEc, vEc, np.log10(error)), ccolors.ANSI_RESET)
if iEc != 0 or vEc != 0:
print(ccolors.ANSI_RED, 'Warning: Solver(s) did not converge', ccolors.ANSI_RESET)
if write:
@@ -143,7 +143,7 @@ class Coupler:
print('{:>5s}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>6s} | {:>5s} {:>5s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'iOut', 'vOut', 'Error'))
print(' --------------------------------------------------------------------')
if couplIter % self.iterPrint == 0:
- print(' {:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>5.0f} {:>5.0f} | {:>6.3f}'.format(couplIter, self.isol.getCl(), self.isol.getTotalDrag(), self.vsol.Cdt, self.vsol.bodies[0].getAvgxtr(0), self.vsol.bodies[0].getAvgxtr(1), iEc, vEc, np.log10(error)))
+ print(' {:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>5.0f} {:>5.0f} | {:>6.3f}'.format(couplIter, self.isol.getCl(), self.isol.getTotalDrag(), self.vsol.getCdt(), self.vsol.bodies[0].getAvgxtr(0), self.vsol.bodies[0].getAvgxtr(1), iEc, vEc, np.log10(error)))
if self.isol.getVerbose() != 0 or self.vsol.verbose != 0:
print('')
couplIter += 1
@@ -152,7 +152,7 @@ class Coupler:
print('')
print('{:>5s}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>8s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
print(' --------------------------------------------------------------------')
- print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter-1, self.isol.getCl(), self.isol.getTotalDrag(), self.vsol.Cdt, self.vsol.bodies[0].getAvgxtr(0), self.vsol.bodies[0].getAvgxtr(1), np.log10(error)), ccolors.ANSI_RESET)
+ print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter-1, self.isol.getCl(), self.isol.getTotalDrag(), self.vsol.getCdt(), self.vsol.bodies[0].getAvgxtr(0), self.vsol.bodies[0].getAvgxtr(1), np.log10(error)), ccolors.ANSI_RESET)
if write:
self.isol.writeCp(sfx='_viscous'+self.filesfx)
for i, body in enumerate(self.vsol.bodies):
diff --git a/blast/blUtils.py b/blast/blUtils.py
index 2b4e24858dd22d491f3423f43bc7c183374470d5..1bfe4d7f05372bb107200173ec2d06fb23d22b49 100644
--- a/blast/blUtils.py
+++ b/blast/blUtils.py
@@ -267,9 +267,9 @@ def save(sections, **kwargs):
region_points.append([nod.pos[0], nod.pos[1], nod.pos[2]])
for var in solution_dict.keys():
- value = getattr(reg, var)[inod]
+ value = float(getattr(reg, var)[inod])
if var == "cf" or var == "cd":
- value *= reg.vt[inod] * reg.vt[inod]
+ value *= float(reg.vt[inod]) * float(reg.vt[inod])
region_solution[var].append(value)
region_points = np.array(region_points)
diff --git a/blast/interfaces/blSolversInterface.py b/blast/interfaces/blSolversInterface.py
index ced73ad27bf6d8725957a3d16a415d69b76852d7..66cec4131f6dd3ef12499bf9337b6df5da1fac71 100644
--- a/blast/interfaces/blSolversInterface.py
+++ b/blast/interfaces/blSolversInterface.py
@@ -122,7 +122,7 @@ class SolversInterface:
float
Total drag coefficient.
"""
- return self.getCd() + self.vSolver.Cdf
+ return self.getCd() + self.vSolver.getCdf()
def getBlowingBoundary(self, couplIter:int):
"""Get blowing boundary condition from the viscous solver.
@@ -361,7 +361,7 @@ class SolversInterface:
# Create node
nodes = blast.std_vector_blNode()
for inod, nod in enumerate(pystruct.getNodes(reg.getName())):
- pos = tbox.Vector3d(nod[xIdx], nod[yIdx], nod[zIdx])
+ pos = blast.blVector3d(nod[xIdx], nod[yIdx], nod[zIdx])
row = int(pystruct.getNodeConnectList(reg.getName())[inod]) if self.cfg['interpolator'] == 'Matching' else inod
nodes.push_back(blast.blNode(inod, pos, row))
# Set nodes
diff --git a/blast/src/CMakeLists.txt b/blast/src/CMakeLists.txt
index cd3f09b324b535c8eaa3e422df8faa6af56dfd51..9d3617aaf25beeedb3888b3607037f1592e986a2 100644
--- a/blast/src/CMakeLists.txt
+++ b/blast/src/CMakeLists.txt
@@ -23,6 +23,11 @@ TARGET_INCLUDE_DIRECTORIES(blast PUBLIC ${PROJECT_SOURCE_DIR}/blast/src
${PROJECT_SOURCE_DIR}/modules/dartflo/ext/amfe/fwk/src
${PROJECT_SOURCE_DIR}/modules/dartflo/dart/src
${PYTHON_INCLUDE_PATH})
+IF(USE_CODI)
+ TARGET_INCLUDE_DIRECTORIES(blast PUBLIC ${CODI_INCLUDE_DIRS})
+ TARGET_COMPILE_DEFINITIONS(blast PUBLIC CODI_EnableEigen)
+ENDIF()
+
TARGET_LINK_LIBRARIES(blast dart fwk tbox)
INSTALL(TARGETS blast DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/blast/src/blBody.cpp b/blast/src/blBody.cpp
index 63c91210a7773c37b95e79e7a2c2c8c48661fdde..d56b5820190e4e01021c1da1e4d4f33cc3719b6f 100644
--- a/blast/src/blBody.cpp
+++ b/blast/src/blBody.cpp
@@ -14,8 +14,8 @@
* limitations under the License.
*/
+#include "blast.h"
#include "blBody.h"
-#include <Eigen/Dense>
using namespace blast;
@@ -27,15 +27,15 @@ void Body::add(std::shared_ptr<Section> section)
sections.push_back(section);
}
-double Body::getAvgxtr(size_t iregion) const
+bldouble Body::getAvgxtr(size_t iregion) const
{
// Return the average transition location of a region.
- double sum = 0.0;
+ bldouble sum = 0.0;
for (auto §ion : sections)
{
sum += section->regions[iregion]->getxoctr();
}
- return sum / static_cast<double>(sections.size());
+ return sum / static_cast<bldouble>(sections.size());
}
void Body::computeDrag()
@@ -48,7 +48,7 @@ void Body::computeDrag()
for (size_t isection = 1; isection < sections.size(); isection++)
{
sections[isection]->computeDrag();
- double dz = sections[isection]->regions[0]->nodes[0]->pos[2] - sections[isection - 1]->regions[0]->nodes[0]->pos[2];
+ bldouble dz = sections[isection]->regions[0]->nodes[0]->pos[2] - sections[isection - 1]->regions[0]->nodes[0]->pos[2];
Cdt += sections[isection]->getCdt() * dz;
Cdf += sections[isection]->getCdf() * dz;
Cdp += sections[isection]->getCdp() * dz;
diff --git a/blast/src/blBody.h b/blast/src/blBody.h
index c868597e8ae60c4de7433e2710fab9a24e9b071f..c029eac4c7c7bd526bc274741cd7310444caa170 100644
--- a/blast/src/blBody.h
+++ b/blast/src/blBody.h
@@ -19,22 +19,22 @@ class BLAST_API Body : public fwk::wSharedObject
{
private:
std::string name; ///< Name of the body.
- double span; ///< Body span. Used for drag computation.
+ bldouble span; ///< Body span. Used for drag computation.
public:
- double Cdt; ///< Total drag coefficient of the body
- double Cdf; ///< Friction drag coefficient of the body
- double Cdp; ///< Pressure drag coefficient of the body
+ bldouble Cdt; ///< Total drag coefficient of the body
+ bldouble Cdf; ///< Friction drag coefficient of the body
+ bldouble Cdp; ///< Pressure drag coefficient of the body
Body(std::string name, double span = 1.0);
~Body() { std::cout << "~blast::Body()" << std::endl; };
std::vector<std::shared_ptr<Section>> sections; ///< Boundary layer sections.
void add(std::shared_ptr<Section> section);
- double getAvgxtr(size_t iregion) const;
+ bldouble getAvgxtr(size_t iregion) const;
void computeDrag();
std::string getName() const { return name; }
- double getSpan() const { return span; }
+ bldouble getSpan() const { return span; }
bool hasWake() const;
};
} // namespace blast
diff --git a/blast/src/blBoundaryLayer.cpp b/blast/src/blBoundaryLayer.cpp
index b30b1e931ddd48dc0e293c190a9ce2ed7a66a2de..9e3ca6e65839236c12bb2e2842f1d1ecfd2325d2 100644
--- a/blast/src/blBoundaryLayer.cpp
+++ b/blast/src/blBoundaryLayer.cpp
@@ -1,3 +1,4 @@
+#include "blast.h"
#include "blBoundaryLayer.h"
#include "blClosures.h"
#include "blNode.h"
@@ -74,9 +75,9 @@ void BoundaryLayer::reset()
*
* @param Re Freestream Reynolds number.
*/
-void BoundaryLayer::setStagnationBC(double Re)
+void BoundaryLayer::setStagnationBC(bldouble Re)
{
- double dv0 = std::abs((vt[1] - vt[0]) / (nodes[1]->xi - nodes[0]->xi));
+ bldouble dv0 = std::abs((vt[1] - vt[0]) / (nodes[1]->xi - nodes[0]->xi));
if (dv0 > 0)
u[0] = sqrt(0.075 / (Re * dv0)); // Theta
else
@@ -94,7 +95,7 @@ void BoundaryLayer::setStagnationBC(double Re)
}
deltaStar[0] = u[0] * u[1];
- std::vector<double> lamParam(8, 0.);
+ std::vector<bldouble> lamParam(8, 0.);
closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name);
Hstar[0] = lamParam[0];
@@ -114,7 +115,7 @@ void BoundaryLayer::setStagnationBC(double Re)
* @param UpperCond Variables at the last point on the upper side.
* @param LowerCond Variables at the last point on the lower side.
*/
-void BoundaryLayer::setWakeBC(double Re, std::vector<std::vector<double>> &teConditions)
+void BoundaryLayer::setWakeBC(bldouble Re, std::vector<std::vector<bldouble>> &teConditions)
{
if (name != "wake")
std::cout << "Warning: Imposing wake boundary condition on " << name
@@ -155,7 +156,7 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<std::vector<double>> &teCon
deltaStar[0] = u[0] * u[1];
// Laminar closures.
- std::vector<double> lamParam(8, 0.);
+ std::vector<bldouble> lamParam(8, 0.);
closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name);
Hstar[0] = lamParam[0];
Hstar2[0] = lamParam[1];
@@ -202,10 +203,10 @@ void BoundaryLayer::computeBlowingVelocity()
}
}
-std::vector<std::vector<double>> BoundaryLayer::evalGradwrtDeltaStar()
+std::vector<std::vector<bldouble>> BoundaryLayer::evalGradwrtDeltaStar()
{
- std::vector<std::vector<double>> gradVe =
- std::vector<std::vector<double>>(getnElms());
+ std::vector<std::vector<bldouble>> gradVe =
+ std::vector<std::vector<bldouble>>(getnElms());
for (auto &vec : gradVe)
vec.resize(getnNodes(), 0.);
@@ -218,10 +219,10 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtDeltaStar()
return gradVe;
}
-std::vector<std::vector<double>> BoundaryLayer::evalGradwrtVt()
+std::vector<std::vector<bldouble>> BoundaryLayer::evalGradwrtVt()
{
- std::vector<std::vector<double>> gradVe =
- std::vector<std::vector<double>>(getnElms());
+ std::vector<std::vector<bldouble>> gradVe =
+ std::vector<std::vector<bldouble>>(getnElms());
for (auto &vec : gradVe)
vec.resize(getnNodes(), 0.);
@@ -234,10 +235,10 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtVt()
return gradVe;
}
-std::vector<std::vector<double>> BoundaryLayer::evalGradwrtRho()
+std::vector<std::vector<bldouble>> BoundaryLayer::evalGradwrtRho()
{
- std::vector<std::vector<double>> gradVe =
- std::vector<std::vector<double>>(getnElms());
+ std::vector<std::vector<bldouble>> gradVe =
+ std::vector<std::vector<bldouble>>(getnElms());
for (auto &vec : gradVe)
vec.resize(getnNodes(), 0.);
@@ -260,25 +261,25 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtRho()
*
* dVe_dx = dVe_dloc * dloc_dx
*
- * @return std::vector<double> Gradient of the velocity with respect to the x
+ * @return std::vector<bldouble> Gradient of the velocity with respect to the x
* coordinate.
*/
-std::vector<std::vector<double>> BoundaryLayer::evalGradwrtX()
+std::vector<std::vector<bldouble>> BoundaryLayer::evalGradwrtX()
{
- std::vector<std::vector<double>> gradVe =
- std::vector<std::vector<double>>(getnElms());
+ std::vector<std::vector<bldouble>> gradVe =
+ std::vector<std::vector<bldouble>>(getnElms());
for (auto &vec : gradVe)
vec.resize(getnNodes(), 0.);
// Compute dVe/dloc
- std::vector<std::vector<double>> gradVe_loc =
- std::vector<std::vector<double>>(getnElms());
+ std::vector<std::vector<bldouble>> gradVe_loc =
+ std::vector<std::vector<bldouble>>(getnElms());
for (auto &vec : gradVe_loc)
vec.resize(getnNodes(), 0.);
for (size_t i = 1; i < getnNodes(); ++i)
{
- double val = 1 / rhoe[i] *
+ bldouble val = 1 / rhoe[i] *
(rhoe[i] * vt[i] * deltaStar[i] -
rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) /
((nodes[i]->xi - nodes[i - 1]->xi) * (nodes[i]->xi - nodes[i - 1]->xi));
@@ -292,8 +293,8 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtX()
// (y_j - y_j-1)^2) = - dloc/dx_j For loc_i, there are two contributions of
// x_j except if j = i. loc_n = sqrt((xn-xn-1)^2) + sqrt((xn-1-xn-2)^2) + ...
// + sqrt((x1-x0)^2)
- std::vector<std::vector<double>> gradloc_x =
- std::vector<std::vector<double>>(getnNodes());
+ std::vector<std::vector<bldouble>> gradloc_x =
+ std::vector<std::vector<bldouble>>(getnNodes());
for (auto &vec : gradloc_x)
vec.resize(getnNodes(), 0.);
@@ -318,22 +319,22 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtX()
return gradVe;
}
-std::vector<std::vector<double>> BoundaryLayer::evalGradwrtY()
+std::vector<std::vector<bldouble>> BoundaryLayer::evalGradwrtY()
{
- std::vector<std::vector<double>> gradVe =
- std::vector<std::vector<double>>(getnElms());
+ std::vector<std::vector<bldouble>> gradVe =
+ std::vector<std::vector<bldouble>>(getnElms());
for (auto &vec : gradVe)
vec.resize(getnNodes(), 0.);
// Compute dVe/dloc
- std::vector<std::vector<double>> gradVe_loc =
- std::vector<std::vector<double>>(getnElms());
+ std::vector<std::vector<bldouble>> gradVe_loc =
+ std::vector<std::vector<bldouble>>(getnElms());
for (auto &vec : gradVe_loc)
vec.resize(getnNodes(), 0.);
for (size_t i = 1; i < getnNodes(); ++i)
{
- double val = 1 / rhoe[i] *
+ bldouble val = 1 / rhoe[i] *
(rhoe[i] * vt[i] * deltaStar[i] -
rhoe[i - 1] * vt[i - 1] * deltaStar[i - 1]) /
((nodes[i]->xi - nodes[i - 1]->xi) * (nodes[i]->xi - nodes[i - 1]->xi));
@@ -347,8 +348,8 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtY()
// (y_j - y_j-1)^2) = - dloc/dy_j For loc_i, there are two contributions of
// y_j except if j = i. loc_n = sqrt((yn-yn-1)^2) + sqrt((yn-1-yn-2)^2) + ...
// + sqrt((y1-y0)^2)
- std::vector<std::vector<double>> gradloc_y =
- std::vector<std::vector<double>>(getnNodes());
+ std::vector<std::vector<bldouble>> gradloc_y =
+ std::vector<std::vector<bldouble>>(getnNodes());
for (auto &vec : gradloc_y)
vec.resize(getnNodes(), 0.);
@@ -373,11 +374,11 @@ std::vector<std::vector<double>> BoundaryLayer::evalGradwrtY()
return gradVe;
}
-std::vector<double> BoundaryLayer::getTeConditions() const
+std::vector<bldouble> BoundaryLayer::getTeConditions() const
{
// Return [theta, H, N, ue, Ct] at the last node.
size_t nodeId = getnNodes() - 1;
- std::vector<double> teCond;
+ std::vector<bldouble> teCond;
for (size_t ivar = 0; ivar < getnVar(); ++ivar)
teCond.push_back(u[nodeId * getnVar() + ivar]);
return teCond;
diff --git a/blast/src/blBoundaryLayer.h b/blast/src/blBoundaryLayer.h
index 04ab90ceea0f7b730433c7c0acfafc8dbb9d24ca..f89f57cc96da4c9e8cd96883ab551e8a99a2c40b 100644
--- a/blast/src/blBoundaryLayer.h
+++ b/blast/src/blBoundaryLayer.h
@@ -22,7 +22,7 @@ class BLAST_API BoundaryLayer : public fwk::wSharedObject
private:
unsigned int const nVar =
5; ///< Number of variables of the partial differential problem.
- double nCrit = 9.0; ///< Critical amplification factor.
+ bldouble nCrit = 9.0; ///< Critical amplification factor.
std::string name; ///< Name of the region.
public:
@@ -33,38 +33,38 @@ public:
std::vector<int> no;
// Boundary layer variables.
- std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct).
- std::vector<double>
+ std::vector<bldouble> u; ///< Solution vector (theta, H, N, ue, Ct).
+ std::vector<bldouble>
Ret; ///< Reynolds number based on the momentum thickness (theta).
- std::vector<double> cd; ///< Local dissipation coefficient.
- std::vector<double> cf; ///< Local friction coefficient.
- std::vector<double>
+ std::vector<bldouble> cd; ///< Local dissipation coefficient.
+ std::vector<bldouble> cf; ///< Local friction coefficient.
+ std::vector<bldouble>
Hstar; ///< Kinetic energy shape parameter (thetaStar/theta).
- std::vector<double>
+ std::vector<bldouble>
Hstar2; ///< Density shape parameter (deltaStarStar/theta).
- std::vector<double> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)).
- std::vector<double>
+ std::vector<bldouble> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)).
+ std::vector<bldouble>
ctEq; ///< Equilibrium shear stress coefficient (turbulent BL).
- std::vector<double> us; ///< Equivalent normalized wall slip velocity.
- std::vector<double> delta; ///< Boundary layer thickness.
- std::vector<double>
+ std::vector<bldouble> us; ///< Equivalent normalized wall slip velocity.
+ std::vector<bldouble> delta; ///< Boundary layer thickness.
+ std::vector<bldouble>
deltaStar; ///< Displacement thickness (int(1-rho*u/rhoe*ue d_eta)).
- std::vector<double> Me; ///< Mach number.
- std::vector<double> vt; ///< Velocity.
- std::vector<double> rhoe; ///< Density.
- std::vector<double> deltaStarExt; ///< Displacement thickness.
- std::vector<double> vtExt; ///< Previous velocity.
+ std::vector<bldouble> Me; ///< Mach number.
+ std::vector<bldouble> vt; ///< Velocity.
+ std::vector<bldouble> rhoe; ///< Density.
+ std::vector<bldouble> deltaStarExt; ///< Displacement thickness.
+ std::vector<bldouble> vtExt; ///< Previous velocity.
- std::vector<double> blowingVelocity; ///< Blowing velocity.
+ std::vector<bldouble> blowingVelocity; ///< Blowing velocity.
// Transition related variables.
size_t transitionNode; ///< Node id of the transition location.
private:
- double xtrF; ///< Forced transition location in the local frame of reference
+ bldouble xtrF; ///< Forced transition location in the local frame of reference
///< (x/c).
- double xtr; ///< Transition location in the local frame of reference (x/c).
+ bldouble xtr; ///< Transition location in the local frame of reference (x/c).
void reset();
@@ -75,34 +75,34 @@ public:
std::string getName() const { return name; }
size_t getnNodes() const { return nodes.size(); }
size_t getnElms() const { return nodes.size() - 1; }
- double getxtr() const { return xtr; }
- double getxoctr() const { return nodes[transitionNode]->xoc; }
- double getxitr() const { return nodes[transitionNode]->xi; }
- double getxtrF() const { return xtrF; }
- double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); }
+ bldouble getxtr() const { return xtr; }
+ bldouble getxoctr() const { return nodes[transitionNode]->xoc; }
+ bldouble getxitr() const { return nodes[transitionNode]->xi; }
+ bldouble getxtrF() const { return xtrF; }
+ bldouble getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); }
size_t getnVar() const { return nVar; }
- double getnCrit() const { return nCrit; }
- std::vector<double> getTeConditions() const;
+ bldouble getnCrit() const { return nCrit; }
+ std::vector<bldouble> getTeConditions() const;
// Setters
void setMesh(std::vector<std::shared_ptr<blNode>> const &_nodes);
void setElementsId(std::vector<int> &nos);
- void setxtrF(double const _xtrF) { xtrF = _xtrF; }
- void setxtr(double const _xtr) { xtr = _xtr; }
- void setnCrit(double const _nCrit) { nCrit = _nCrit; }
+ void setxtrF(bldouble const _xtrF) { xtrF = _xtrF; }
+ void setxtr(bldouble const _xtr) { xtr = _xtr; }
+ void setnCrit(bldouble const _nCrit) { nCrit = _nCrit; }
// Boundary conditions.
- void setStagnationBC(double Re);
- void setWakeBC(double Re, std::vector<std::vector<double>> &teConditions);
+ void setStagnationBC(bldouble Re);
+ void setWakeBC(bldouble Re, std::vector<std::vector<bldouble>> &teConditions);
// Others
void printSolution(size_t inod) const;
void computeBlowingVelocity();
- std::vector<std::vector<double>> evalGradwrtRho();
- std::vector<std::vector<double>> evalGradwrtVt();
- std::vector<std::vector<double>> evalGradwrtDeltaStar();
- std::vector<std::vector<double>> evalGradwrtX();
- std::vector<std::vector<double>> evalGradwrtY();
+ std::vector<std::vector<bldouble>> evalGradwrtRho();
+ std::vector<std::vector<bldouble>> evalGradwrtVt();
+ std::vector<std::vector<bldouble>> evalGradwrtDeltaStar();
+ std::vector<std::vector<bldouble>> evalGradwrtX();
+ std::vector<std::vector<bldouble>> evalGradwrtY();
};
} // namespace blast
#endif // BLBOUNDARYLAYER_H
diff --git a/blast/src/blClosures.cpp b/blast/src/blClosures.cpp
index e2a4f2ab098bb2fc8a9818de261e20183a832cd7..249b15391a194d87dec0742ddf3cb72f97e14787 100644
--- a/blast/src/blClosures.cpp
+++ b/blast/src/blClosures.cpp
@@ -1,3 +1,4 @@
+#include "blast.h"
#include "blClosures.h"
#include <cmath>
#include <iostream>
@@ -14,27 +15,27 @@ using namespace blast;
* @param Me Local Mach number.
* @param name Name of the region.
*/
-void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
- double H, double Ret, const double Me,
+void Closures::laminarClosures(std::vector<bldouble> &closureVars, bldouble theta,
+ bldouble H, bldouble Ret, const bldouble Me,
const std::string name)
{
H = std::max(H, 1.00005);
// Kinematic shape factor H*.
- double Hk = (H - 0.29 * Me * Me) / (1.0 + 0.113 * Me * Me);
+ bldouble Hk = (H - 0.29 * Me * Me) / (1.0 + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
Hk = std::max(Hk, 1.05000);
// Density shape parameter.
- double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me;
+ bldouble Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me;
// Boundary layer thickness.
- double delta = std::min((3.15 + H + (1.72 / (Hk - 1.0))) * theta, 12.0 * theta);
+ bldouble delta = std::min((3.15 + H + (1.72 / (Hk - 1.0))) * theta, 12.0 * theta);
- double Hstar = 0.0;
- double Hks = Hk - 4.35;
+ bldouble Hstar = 0.0;
+ bldouble Hks = Hk - 4.35;
if (Hk <= 4.35)
Hstar = 0.0111 * Hks * Hks / (Hk + 1.0) -
0.0278 * Hks * Hks * Hks / (Hk + 1.0) + 1.528 -
@@ -46,8 +47,8 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
Hstar = (Hstar + 0.028 * Me * Me) / (1.0 + 0.014 * Me * Me);
// Friction coefficient.
- double tmp = 0.0;
- double cf = 0.0;
+ bldouble tmp = 0.0;
+ bldouble cf = 0.0;
if (Hk < 5.5)
{
tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0);
@@ -60,13 +61,13 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
}
// Dissipation coefficient.
- double Cd = 0.0;
+ bldouble Cd = 0.0;
if (Hk < 4.0)
Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2.0 * Ret));
else
{
- double HkCd = (Hk - 4.0) * (Hk - 4.0);
- double denCd = 1.0 + 0.02 * HkCd;
+ bldouble HkCd = (Hk - 4.0) * (Hk - 4.0);
+ bldouble denCd = 1.0 + 0.02 * HkCd;
Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2.0 * Ret));
}
@@ -78,8 +79,8 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
cf = 0.0;
}
- double us = 0.;
- double ctEq = 0.;
+ bldouble us = 0.;
+ bldouble ctEq = 0.;
closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
}
@@ -95,25 +96,25 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
* @param Me Local Mach number.
* @param name Name of the region.
*/
-void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
- double H, double Ct, double Ret,
- const double Me, const std::string name)
+void Closures::turbulentClosures(std::vector<bldouble> &closureVars, bldouble theta,
+ bldouble H, bldouble Ct, bldouble Ret,
+ const bldouble Me, const std::string name)
{
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
// Ct = std::max(std::min(Ct, 0.30), 0.0000001);
- double Hk = (H - 0.29 * Me * Me) / (1. + 0.113 * Me * Me);
+ bldouble Hk = (H - 0.29 * Me * Me) / (1. + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
Hk = std::max(Hk, 1.05000);
- double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
- double gamma = 1.4 - 1.;
- double Fc = std::sqrt(1. + 0.5 * gamma * Me * Me);
+ bldouble Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
+ bldouble gamma = 1.4 - 1.;
+ bldouble Fc = std::sqrt(1. + 0.5 * gamma * Me * Me);
- double H0 = 0.;
+ bldouble H0 = 0.;
if (Ret <= 400.)
H0 = 4.0;
else
@@ -121,7 +122,7 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
if (Ret <= 200.)
Ret = 200.;
- double Hstar = 0.;
+ bldouble Hstar = 0.;
if (Hk <= H0)
Hstar =
((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
@@ -138,32 +139,32 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
// Whitfield's minor additional compressibility correction.
Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me);
- double logRt = std::log(Ret / Fc);
+ bldouble logRt = std::log(Ret / Fc);
logRt = std::max(logRt, 3.0);
- double arg = std::max(-1.33 * Hk, -20.0);
+ bldouble arg = std::max(-1.33 * Hk, -20.0);
// Equivalent normalized wall slip velocity.
- double us = (Hstar / 2.) * (1. - 4. * (Hk - 1.) / (3. * H));
+ bldouble us = (Hstar / 2.) * (1. - 4. * (Hk - 1.) / (3. * H));
// Boundary layer thickness.
- double delta = std::min((3.15 + H + (1.72 / (Hk - 1.))) * theta, 12. * theta);
+ bldouble delta = std::min((3.15 + H + (1.72 / (Hk - 1.))) * theta, 12. * theta);
- double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
- double Hkc = 0.;
- double Cdw = 0.;
- double Cdd = 0.;
- double Cdl = 0.;
+ bldouble Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
+ bldouble Hkc = 0.;
+ bldouble Cdw = 0.;
+ bldouble Cdd = 0.;
+ bldouble Cdl = 0.;
- double Hmin = 0.;
- double Fl = 0.;
- double Dfac = 0.;
+ bldouble Hmin = 0.;
+ bldouble Fl = 0.;
+ bldouble Dfac = 0.;
- double cf = 0.;
- double Cd = 0.;
+ bldouble cf = 0.;
+ bldouble Cd = 0.;
- double n = 1.0;
+ bldouble n = 1.0;
- double ctEq = 0.;
+ bldouble ctEq = 0.;
if (name == "wake")
{
@@ -221,20 +222,20 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
* @param Me Local Mach number.
* @param name Name of the region.
*/
-void Closures::turbulentClosures(double &closureVars, double theta, double H,
- double Ct, double Ret, const double Me,
+void Closures::turbulentClosures(bldouble &closureVars, bldouble theta, bldouble H,
+ bldouble Ct, bldouble Ret, const bldouble Me,
const std::string name)
{
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
- double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
+ bldouble Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
Hk = std::max(Hk, 1.05000);
- double H0 = 0.0;
+ bldouble H0 = 0.0;
if (Ret <= 400.)
H0 = 4.0;
else
@@ -242,7 +243,7 @@ void Closures::turbulentClosures(double &closureVars, double theta, double H,
if (Ret <= 200.0)
Ret = 200.0;
- double Hstar = 0.0;
+ bldouble Hstar = 0.0;
if (Hk <= H0)
Hstar =
((0.5 - 4.0 / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
@@ -260,12 +261,12 @@ void Closures::turbulentClosures(double &closureVars, double theta, double H,
Hstar = (Hstar + 0.028 * Me * Me) / (1.0 + 0.014 * Me * Me);
// Equivalent normalized wall slip velocity.
- double us = (Hstar / 2.0) * (1.0 - 4 * (Hk - 1.0) / (3.0 * H));
+ bldouble us = (Hstar / 2.0) * (1.0 - 4 * (Hk - 1.0) / (3.0 * H));
- double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
+ bldouble Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
- double Hkc = 0.0;
- double ctEq = 0.0;
+ bldouble Hkc = 0.0;
+ bldouble ctEq = 0.0;
if (name == "wake")
{
diff --git a/blast/src/blClosures.h b/blast/src/blClosures.h
index 42d48ecfcc1d61bf3267c3240ca156294b360a7a..7fd6f4cc3e70bdb1d067810f2dd641f16df3f83d 100644
--- a/blast/src/blClosures.h
+++ b/blast/src/blClosures.h
@@ -16,14 +16,14 @@ class BLAST_API Closures
{
public:
- static void laminarClosures(std::vector<double> &closureVars, double theta,
- double H, double Ret, const double Me,
+ static void laminarClosures(std::vector<bldouble> &closureVars, bldouble theta,
+ bldouble H, bldouble Ret, const bldouble Me,
const std::string name);
- static void turbulentClosures(std::vector<double> &closureVars, double theta,
- double H, double Ct, double Ret, double Me,
+ static void turbulentClosures(std::vector<bldouble> &closureVars, bldouble theta,
+ bldouble H, bldouble Ct, bldouble Ret, bldouble Me,
const std::string name);
- static void turbulentClosures(double &closureVars, double theta, double H,
- double Ct, double Ret, const double Me,
+ static void turbulentClosures(bldouble &closureVars, bldouble theta, bldouble H,
+ bldouble Ct, bldouble Ret, const bldouble Me,
const std::string name);
};
} // namespace blast
diff --git a/blast/src/blCoupledAdjoint.cpp b/blast/src/blCoupledAdjoint.cpp
index 55421a7c5da2607c924db5911344f77d9a5602a8..1e8867fbadb54dfa0fc7388d3c34cf16b8eb0b54 100644
--- a/blast/src/blCoupledAdjoint.cpp
+++ b/blast/src/blCoupledAdjoint.cpp
@@ -576,97 +576,97 @@ void CoupledAdjoint::gradientswrtMachNumber()
//**************************************************************************//
dRM_M.setIdentity();
- //**************************************************************************//
- // dRdStar_M //
- //--------------------------------------------------------------------------//
- // Central finite-difference //
- //**************************************************************************//
- std::deque<Eigen::Triplet<double>> T;
- double delta = 1e-6;
- size_t i = 0;
- double saveM = 0.;
- std::vector<Eigen::VectorXd> ddStar(2,
- Eigen::VectorXd::Zero(dRdStar_M.cols()));
- std::vector<double> dCdf(2, 0.);
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- for (size_t j = 0; j < reg->getnNodes(); ++j)
- {
- saveM = reg->Me[j];
-
- reg->Me[j] = saveM + delta;
- ddStar[0] = this->runViscous();
- dCdf[0] = vsol->Cdf;
- reg->Me[j] = saveM - delta;
- ddStar[1] = this->runViscous();
- dCdf[1] = vsol->Cdf;
-
- reg->Me[j] = saveM;
-
- for (int k = 0; k < ddStar[0].size(); ++k)
- T.push_back(Eigen::Triplet<double>(
- k, i, -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
- // Collect Cd gradients
- dCd_M(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
- ++i;
- }
- }
- }
- dRdStar_M.setFromTriplets(T.begin(), T.end());
- dRdStar_M.prune(0.);
- }
+// //**************************************************************************//
+// // dRdStar_M //
+// //--------------------------------------------------------------------------//
+// // Central finite-difference //
+// //**************************************************************************//
+// std::deque<Eigen::Triplet<double>> T;
+// double delta = 1e-6;
+// size_t i = 0;
+// double saveM = 0.;
+// std::vector<Eigen::VectorXd> ddStar(2,
+// Eigen::VectorXd::Zero(dRdStar_M.cols()));
+// std::vector<double> dCdf(2, 0.);
+// for (auto &body : vsol->bodies)
+// {
+// for (auto §ion : body->sections)
+// {
+// for (auto ® : section->regions)
+// {
+// for (size_t j = 0; j < reg->getnNodes(); ++j)
+// {
+// saveM = reg->Me[j];
+
+// reg->Me[j] = saveM + delta;
+// ddStar[0] = this->runViscous();
+// dCdf[0] = vsol->Cdf;
+// reg->Me[j] = saveM - delta;
+// ddStar[1] = this->runViscous();
+// dCdf[1] = vsol->Cdf;
+
+// reg->Me[j] = saveM;
+
+// for (int k = 0; k < ddStar[0].size(); ++k)
+// T.push_back(Eigen::Triplet<double>(
+// k, i, -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
+// // Collect Cd gradients
+// dCd_M(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
+// ++i;
+// }
+// }
+// }
+// dRdStar_M.setFromTriplets(T.begin(), T.end());
+// dRdStar_M.prune(0.);
+// }
}
-/**
- * @brief Compute all the gradients wrt the density on the surface
- * Non-zero are: dRrho_rho, dRblw_rho
- * Zeros are: dRphi_rho, dRM_rho, dRv_rho, dRdStar_rho, dRx_rho
- */
+// /**
+// * @brief Compute all the gradients wrt the density on the surface
+// * Non-zero are: dRrho_rho, dRblw_rho
+// * Zeros are: dRphi_rho, dRM_rho, dRv_rho, dRdStar_rho, dRx_rho
+// */
void CoupledAdjoint::gradientswrtDensity()
{
- size_t nNs = 0;
- for (auto bBC : isol->pbl->bBCs)
- nNs += bBC->nodes.size(); // Number of surface nodes
- nNs += 1; // Duplicate stagnation point
- size_t nEs = 0;
- for (auto bBC : isol->pbl->bBCs)
- nEs += bBC->uE.size(); // Number of blowing elements
-
- //**************************************************************************//
- // dRrho_rho //
- //**************************************************************************//
- dRrho_rho.setIdentity();
-
- //**************************************************************************//
- // dRblw_rho //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- std::deque<Eigen::Triplet<double>> T_blw;
- int offSetElms = 0;
- int offSetNodes = 0;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- std::vector<std::vector<double>> grad = reg->evalGradwrtRho();
- for (size_t i = 0; i < grad.size(); ++i)
- for (size_t j = 0; j < grad[i].size(); ++j)
- T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
- -grad[i][j]));
- offSetElms += reg->getnElms();
- offSetNodes += reg->getnNodes();
- }
- }
- }
- dRblw_rho.setFromTriplets(T_blw.begin(), T_blw.end());
- dRblw_rho.prune(0.);
+// size_t nNs = 0;
+// for (auto bBC : isol->pbl->bBCs)
+// nNs += bBC->nodes.size(); // Number of surface nodes
+// nNs += 1; // Duplicate stagnation point
+// size_t nEs = 0;
+// for (auto bBC : isol->pbl->bBCs)
+// nEs += bBC->uE.size(); // Number of blowing elements
+
+// //**************************************************************************//
+// // dRrho_rho //
+// //**************************************************************************//
+// dRrho_rho.setIdentity();
+
+// //**************************************************************************//
+// // dRblw_rho //
+// //--------------------------------------------------------------------------//
+// // Analytical //
+// //**************************************************************************//
+// std::deque<Eigen::Triplet<double>> T_blw;
+// int offSetElms = 0;
+// int offSetNodes = 0;
+// for (auto &body : vsol->bodies)
+// {
+// for (auto §ion : body->sections)
+// {
+// for (auto ® : section->regions)
+// {
+// std::vector<std::vector<double>> grad = reg->evalGradwrtRho();
+// for (size_t i = 0; i < grad.size(); ++i)
+// for (size_t j = 0; j < grad[i].size(); ++j)
+// T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
+// -grad[i][j]));
+// offSetElms += reg->getnElms();
+// offSetNodes += reg->getnNodes();
+// }
+// }
+// }
+// dRblw_rho.setFromTriplets(T_blw.begin(), T_blw.end());
+// dRblw_rho.prune(0.);
}
/**
@@ -676,168 +676,168 @@ void CoupledAdjoint::gradientswrtDensity()
*/
void CoupledAdjoint::gradientswrtVelocity()
{
- size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
- size_t nNs = 0;
- for (auto bBC : isol->pbl->bBCs)
- nNs += bBC->nodes.size(); // Number of surface nodes
- nNs += 1; // Duplicate stagnation point
- size_t nEs = 0;
- for (auto bBC : isol->pbl->bBCs)
- nEs += bBC->uE.size(); // Number of blowing elements
-
- //**************************************************************************//
- // dRv_v //
- //**************************************************************************//
- dRv_v.setIdentity();
-
- //**** Get velocity****/
- Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nNs, nDim);
- size_t i = 0;
- int iReg = -1;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- if (reg->getName() == "wake")
- iReg = 1;
- else if (reg->getName() == "upper" || reg->getName() == "lower")
- iReg = 0;
- else
- throw std::runtime_error("Unknown section name");
- for (size_t j = 0; j < reg->getnNodes(); ++j)
- {
- for (size_t idim = 0; idim < nDim; ++idim)
- v(i, idim) =
- isol->U[isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row](idim);
- ++i;
- }
- }
- }
- }
-
- //**************************************************************************//
- // dvt_v //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim);
- std::deque<Eigen::Triplet<double>> T_vt;
-
- i = 0;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- for (size_t j = 0; j < reg->getnNodes(); ++j)
- {
- for (size_t idim = 0; idim < nDim; ++idim)
- {
- T_vt.push_back(Eigen::Triplet<double>(
- i, i * nDim + idim,
- 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim)));
- }
- ++i;
- }
- }
- }
- }
- dVt_v.setFromTriplets(T_vt.begin(), T_vt.end());
- dVt_v.prune(0.);
-
- //**************************************************************************//
- // dRdStar_vt //
- //--------------------------------------------------------------------------//
- // Central finite-difference //
- //**************************************************************************//
- std::deque<Eigen::Triplet<double>> T;
- double delta = 1e-6;
- i = 0;
- double saveM = 0.;
- std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
- std::vector<double> dCdf(2, 0.);
- Eigen::VectorXd dCd_vt = Eigen::VectorXd::Zero(dRM_M.cols());
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- for (size_t j = 0; j < reg->getnNodes(); ++j)
- {
- saveM = reg->vt[j];
-
- reg->vt[j] = saveM + delta;
- dvt[0] = this->runViscous();
- dCdf[0] = vsol->Cdf;
- reg->vt[j] = saveM - delta;
- dvt[1] = this->runViscous();
- dCdf[1] = vsol->Cdf;
-
- reg->vt[j] = saveM;
-
- for (int k = 0; k < dvt[0].size(); ++k)
- T.push_back(Eigen::Triplet<double>(
- k, i, -(dvt[0][k] - dvt[1][k]) / (2. * delta)));
- dCd_vt(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
- ++i;
- }
- }
- }
- }
- Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(),
- dRdStar_M.cols());
- dRdStar_vt.setFromTriplets(T.begin(), T.end());
- dRdStar_vt.prune(0.);
-
- //**************************************************************************//
- // dRdStar_v //
- //**************************************************************************//
- dRdStar_v = dRdStar_vt * dVt_v;
- dRdStar_v.prune(0.);
-
- //**************************************************************************//
- // dCdf_v //
- //**************************************************************************//
- dCd_v =
- dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim]
-
- //**************************************************************************//
- // dRblw_vt //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- std::deque<Eigen::Triplet<double>> T_blw;
- int offSetElms = 0;
- int offSetNodes = 0;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- std::vector<std::vector<double>> grad = reg->evalGradwrtVt();
- for (size_t i = 0; i < grad.size(); ++i)
- for (size_t j = 0; j < grad[i].size(); ++j)
- T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
- -grad[i][j]));
- offSetElms += reg->getnElms();
- offSetNodes += reg->getnNodes();
- }
- }
- }
- Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt(nEs, nNs);
- dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end());
- dRblw_vt.prune(0.);
-
- //**************************************************************************//
- // dRblw_v //
- //**************************************************************************//
- dRblw_v = dRblw_vt * dVt_v;
- dRblw_v.prune(0.);
+ // size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
+ // size_t nNs = 0;
+ // for (auto bBC : isol->pbl->bBCs)
+ // nNs += bBC->nodes.size(); // Number of surface nodes
+ // nNs += 1; // Duplicate stagnation point
+ // size_t nEs = 0;
+ // for (auto bBC : isol->pbl->bBCs)
+ // nEs += bBC->uE.size(); // Number of blowing elements
+
+ // //**************************************************************************//
+ // // dRv_v //
+ // //**************************************************************************//
+ // dRv_v.setIdentity();
+
+ // //**** Get velocity****/
+ // Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nNs, nDim);
+ // size_t i = 0;
+ // int iReg = -1;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // if (reg->getName() == "wake")
+ // iReg = 1;
+ // else if (reg->getName() == "upper" || reg->getName() == "lower")
+ // iReg = 0;
+ // else
+ // throw std::runtime_error("Unknown section name");
+ // for (size_t j = 0; j < reg->getnNodes(); ++j)
+ // {
+ // for (size_t idim = 0; idim < nDim; ++idim)
+ // v(i, idim) =
+ // isol->U[isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row](idim);
+ // ++i;
+ // }
+ // }
+ // }
+ // }
+
+ // //**************************************************************************//
+ // // dvt_v //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim);
+ // std::deque<Eigen::Triplet<double>> T_vt;
+
+ // i = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // for (size_t j = 0; j < reg->getnNodes(); ++j)
+ // {
+ // for (size_t idim = 0; idim < nDim; ++idim)
+ // {
+ // T_vt.push_back(Eigen::Triplet<double>(
+ // i, i * nDim + idim,
+ // 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim)));
+ // }
+ // ++i;
+ // }
+ // }
+ // }
+ // }
+ // dVt_v.setFromTriplets(T_vt.begin(), T_vt.end());
+ // dVt_v.prune(0.);
+
+ // //**************************************************************************//
+ // // dRdStar_vt //
+ // //--------------------------------------------------------------------------//
+ // // Central finite-difference //
+ // //**************************************************************************//
+ // std::deque<Eigen::Triplet<double>> T;
+ // double delta = 1e-6;
+ // i = 0;
+ // double saveM = 0.;
+ // std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ // std::vector<double> dCdf(2, 0.);
+ // Eigen::VectorXd dCd_vt = Eigen::VectorXd::Zero(dRM_M.cols());
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // for (size_t j = 0; j < reg->getnNodes(); ++j)
+ // {
+ // saveM = reg->vt[j];
+
+ // reg->vt[j] = saveM + delta;
+ // dvt[0] = this->runViscous();
+ // dCdf[0] = vsol->Cdf;
+ // reg->vt[j] = saveM - delta;
+ // dvt[1] = this->runViscous();
+ // dCdf[1] = vsol->Cdf;
+
+ // reg->vt[j] = saveM;
+
+ // for (int k = 0; k < dvt[0].size(); ++k)
+ // T.push_back(Eigen::Triplet<double>(
+ // k, i, -(dvt[0][k] - dvt[1][k]) / (2. * delta)));
+ // dCd_vt(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
+ // ++i;
+ // }
+ // }
+ // }
+ // }
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(),
+ // dRdStar_M.cols());
+ // dRdStar_vt.setFromTriplets(T.begin(), T.end());
+ // dRdStar_vt.prune(0.);
+
+ // //**************************************************************************//
+ // // dRdStar_v //
+ // //**************************************************************************//
+ // dRdStar_v = dRdStar_vt * dVt_v;
+ // dRdStar_v.prune(0.);
+
+ // //**************************************************************************//
+ // // dCdf_v //
+ // //**************************************************************************//
+ // dCd_v =
+ // dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim]
+
+ // //**************************************************************************//
+ // // dRblw_vt //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // std::deque<Eigen::Triplet<double>> T_blw;
+ // int offSetElms = 0;
+ // int offSetNodes = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // std::vector<std::vector<double>> grad = reg->evalGradwrtVt();
+ // for (size_t i = 0; i < grad.size(); ++i)
+ // for (size_t j = 0; j < grad[i].size(); ++j)
+ // T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
+ // -grad[i][j]));
+ // offSetElms += reg->getnElms();
+ // offSetNodes += reg->getnNodes();
+ // }
+ // }
+ // }
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt(nEs, nNs);
+ // dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end());
+ // dRblw_vt.prune(0.);
+
+ // //**************************************************************************//
+ // // dRblw_v //
+ // //**************************************************************************//
+ // dRblw_v = dRblw_vt * dVt_v;
+ // dRblw_v.prune(0.);
}
/**
@@ -852,77 +852,77 @@ void CoupledAdjoint::gradientswrtDeltaStar()
//--------------------------------------------------------------------------//
// Central finite-difference //
//**************************************************************************//
- std::deque<Eigen::Triplet<double>> T;
- double delta = 1e-6;
- double saveDStar = 0.;
- std::vector<Eigen::VectorXd> ddstar(
- 2, Eigen::VectorXd::Zero(dRdStar_dStar.cols()));
- std::vector<double> dCdf(2, 0.);
- size_t i = 0;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- for (size_t j = 0; j < reg->getnNodes(); ++j)
- {
- saveDStar = reg->deltaStarExt[j];
-
- reg->deltaStarExt[j] = saveDStar + delta;
- ddstar[0] = this->runViscous();
- dCdf[0] = vsol->Cdf;
- reg->deltaStarExt[j] = saveDStar - delta;
- ddstar[1] = this->runViscous();
- dCdf[1] = vsol->Cdf;
-
- reg->deltaStarExt[j] = saveDStar;
-
- for (int k = 0; k < ddstar[0].size(); ++k)
- T.push_back(Eigen::Triplet<double>(
- k, i, (ddstar[0][k] - ddstar[1][k]) / (2. * delta)));
- dCd_dStar(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
- ++i;
- }
- }
- }
- }
- Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum(
- dRdStar_dStar.rows(), dRdStar_dStar.cols());
- // RdStar = dStar - f(dStar) = 0
- // dRdStar_dStar = 1 - df_dStar (minus is 3 lines below and not in the
- // triplets)
- dRdStar_dStar_dum.setFromTriplets(T.begin(), T.end());
- dRdStar_dStar.setIdentity();
- dRdStar_dStar -= dRdStar_dStar_dum;
- dRdStar_dStar.prune(0.);
-
- //**************************************************************************//
- // dRblw_dStar //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- std::deque<Eigen::Triplet<double>> T_blw;
- int offSetElms = 0;
- int offSetNodes = 0;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- std::vector<std::vector<double>> grad = reg->evalGradwrtDeltaStar();
- for (size_t i = 0; i < grad.size(); ++i)
- for (size_t j = 0; j < grad[i].size(); ++j)
- T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
- -grad[i][j]));
- offSetElms += reg->getnElms();
- offSetNodes += reg->getnNodes();
- }
- }
- }
- dRblw_dStar.setFromTriplets(T_blw.begin(), T_blw.end());
- dRblw_dStar.prune(0.);
+ // std::deque<Eigen::Triplet<double>> T;
+ // double delta = 1e-6;
+ // double saveDStar = 0.;
+ // std::vector<Eigen::VectorXd> ddstar(
+ // 2, Eigen::VectorXd::Zero(dRdStar_dStar.cols()));
+ // std::vector<double> dCdf(2, 0.);
+ // size_t i = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // for (size_t j = 0; j < reg->getnNodes(); ++j)
+ // {
+ // saveDStar = reg->deltaStarExt[j];
+
+ // reg->deltaStarExt[j] = saveDStar + delta;
+ // ddstar[0] = this->runViscous();
+ // dCdf[0] = vsol->Cdf;
+ // reg->deltaStarExt[j] = saveDStar - delta;
+ // ddstar[1] = this->runViscous();
+ // dCdf[1] = vsol->Cdf;
+
+ // reg->deltaStarExt[j] = saveDStar;
+
+ // for (int k = 0; k < ddstar[0].size(); ++k)
+ // T.push_back(Eigen::Triplet<double>(
+ // k, i, (ddstar[0][k] - ddstar[1][k]) / (2. * delta)));
+ // dCd_dStar(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
+ // ++i;
+ // }
+ // }
+ // }
+ // }
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum(
+ // dRdStar_dStar.rows(), dRdStar_dStar.cols());
+ // // RdStar = dStar - f(dStar) = 0
+ // // dRdStar_dStar = 1 - df_dStar (minus is 3 lines below and not in the
+ // // triplets)
+ // dRdStar_dStar_dum.setFromTriplets(T.begin(), T.end());
+ // dRdStar_dStar.setIdentity();
+ // dRdStar_dStar -= dRdStar_dStar_dum;
+ // dRdStar_dStar.prune(0.);
+
+ // //**************************************************************************//
+ // // dRblw_dStar //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // std::deque<Eigen::Triplet<double>> T_blw;
+ // int offSetElms = 0;
+ // int offSetNodes = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // std::vector<std::vector<double>> grad = reg->evalGradwrtDeltaStar();
+ // for (size_t i = 0; i < grad.size(); ++i)
+ // for (size_t j = 0; j < grad[i].size(); ++j)
+ // T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
+ // -grad[i][j]));
+ // offSetElms += reg->getnElms();
+ // offSetNodes += reg->getnNodes();
+ // }
+ // }
+ // }
+ // dRblw_dStar.setFromTriplets(T_blw.begin(), T_blw.end());
+ // dRblw_dStar.prune(0.);
}
/**
@@ -937,44 +937,44 @@ void CoupledAdjoint::gradientswrtBlowingVelocity()
//--------------------------------------------------------------------------//
// Analytical //
//**************************************************************************//
- std::deque<Eigen::Triplet<double>> T;
-
- size_t jj = 0;
- int iReg = 0;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- if (reg->getName() == "wake")
- iReg = 1;
- else if (reg->getName() == "upper" || reg->getName() == "lower")
- iReg = 0;
- else
- throw std::runtime_error("Unknown section name");
- for (size_t iElm = 0; iElm < reg->getnElms(); ++iElm)
- {
- auto blowingElement = isol->pbl->bBCs[iReg]->uE[reg->no[iElm]].first;
- Eigen::VectorXd be =
- dart::BlowingResidual::buildGradientBlowing(*blowingElement);
- for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++)
- {
- tbox::Node *nodi = blowingElement->nodes[ii];
- T.push_back(
- Eigen::Triplet<double>(isol->getRow(nodi->row), jj, be(ii)));
- }
- ++jj;
- }
- }
- }
- }
- dRphi_blw.setFromTriplets(T.begin(), T.end());
- dRphi_blw.prune(0.);
- dRphi_blw.makeCompressed();
-
- /**** dRblw_blw ****/
- dRblw_blw.setIdentity();
+ // std::deque<Eigen::Triplet<double>> T;
+
+ // size_t jj = 0;
+ // int iReg = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // if (reg->getName() == "wake")
+ // iReg = 1;
+ // else if (reg->getName() == "upper" || reg->getName() == "lower")
+ // iReg = 0;
+ // else
+ // throw std::runtime_error("Unknown section name");
+ // for (size_t iElm = 0; iElm < reg->getnElms(); ++iElm)
+ // {
+ // auto blowingElement = isol->pbl->bBCs[iReg]->uE[reg->no[iElm]].first;
+ // Eigen::VectorXd be =
+ // dart::BlowingResidual::buildGradientBlowing(*blowingElement);
+ // for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++)
+ // {
+ // tbox::Node *nodi = blowingElement->nodes[ii];
+ // T.push_back(
+ // Eigen::Triplet<double>(isol->getRow(nodi->row), jj, be(ii)));
+ // }
+ // ++jj;
+ // }
+ // }
+ // }
+ // }
+ // dRphi_blw.setFromTriplets(T.begin(), T.end());
+ // dRphi_blw.prune(0.);
+ // dRphi_blw.makeCompressed();
+
+ // /**** dRblw_blw ****/
+ // dRblw_blw.setIdentity();
}
/**
@@ -989,16 +989,16 @@ void CoupledAdjoint::gradientwrtAoA()
//--------------------------------------------------------------------------//
// Analytical //
//**************************************************************************//
- std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa();
- dCl_AoA = dCx_AoA[0];
- dCd_AoA = dCx_AoA[1];
-
- //**************************************************************************//
- // dRphi_Aoa //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- dRphi_AoA = iadjoint->getRu_A();
+ // std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa();
+ // dCl_AoA = dCx_AoA[0];
+ // dCd_AoA = dCx_AoA[1];
+
+ // //**************************************************************************//
+ // // dRphi_Aoa //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // dRphi_AoA = iadjoint->getRu_A();
}
/**
@@ -1009,228 +1009,228 @@ void CoupledAdjoint::gradientwrtAoA()
void CoupledAdjoint::gradientwrtMesh()
{
- int nDim = isol->pbl->nDim;
-
- //**************************************************************************//
- // dCl_x, dCd_x //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- std::vector<std::vector<double>> dCx_x =
- iadjoint->computeGradientCoefficientsMesh();
- dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size());
- dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
-
- //**************************************************************************//
- // dRphi_x //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- dRphi_x = iadjoint->getRu_X();
- //**************************************************************************//
- // dRx_x //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- dRx_x = iadjoint->getRx_X();
-
- //**************************************************************************//
- // dRdStar_x, dCdf_dx //
- //--------------------------------------------------------------------------//
- // Central finite-difference //
- //--------------------------------------------------------------------------//
- // We don't do dRdStar_x = dRdStar_loc * dloc_dx because Cdf = f(xi(x), x) //
- // and we would have to compute dCdf/dx and dCdf/dloc. Same for deltaStar. //
- // It is more computationally expansive to compute //
- // ddeltaStar/dx, ddeltaStar/dloc than ddeltaStar_dx, ddeltaStar_dy //
- //**************************************************************************//
- std::deque<Eigen::Triplet<double>> T_dStar_x;
- double delta = 1e-8;
- int iReg = 0;
- std::vector<Eigen::VectorXd> ddStar(2,
- Eigen::VectorXd::Zero(dRdStar_M.cols()));
- std::vector<double> dCdf(2, 0.);
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- if (reg->getName() == "wake")
- iReg = 1;
- else if (reg->getName() == "upper" || reg->getName() == "lower")
- iReg = 0;
- else
- throw std::runtime_error("Unknown section name");
- for (size_t j = 0; j < reg->getnNodes(); ++j)
- {
- int rowj = isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row;
- // x
- double saveX = reg->nodes[j]->pos[0];
-
- reg->nodes[j]->pos[0] = saveX + delta;
- section->computeLocalCoordinates(reg);
- for (auto &nod : reg->nodes)
- nod->xiExt = nod->xi;
- ddStar[0] = this->runViscous();
- dCdf[0] = vsol->Cdf;
- reg->nodes[j]->pos[0] = saveX - delta;
- section->computeLocalCoordinates(reg);
- for (auto &nod : reg->nodes)
- nod->xiExt = nod->xi;
- ddStar[1] = this->runViscous();
- dCdf[1] = vsol->Cdf;
-
- reg->nodes[j]->pos[0] = saveX;
- section->computeLocalCoordinates(reg);
- for (auto &nod : reg->nodes)
- nod->xiExt = nod->xi;
-
- for (int k = 0; k < ddStar[0].size(); ++k)
- T_dStar_x.push_back(Eigen::Triplet<double>(
- k, rowj * isol->pbl->nDim + 0,
- -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
- dCdf_x(rowj * isol->pbl->nDim + 0) += (dCdf[0] - dCdf[1]) / (2. * delta);
- // y
- double saveY = reg->nodes[j]->pos[1];
-
- reg->nodes[j]->pos[1] = saveY + delta;
- section->computeLocalCoordinates(reg);
- for (auto &nod : reg->nodes)
- nod->xiExt = nod->xi;
- ddStar[0] = this->runViscous();
- dCdf[0] = vsol->Cdf;
- reg->nodes[j]->pos[1] = saveY - delta;
- section->computeLocalCoordinates(reg);
- for (auto &nod : reg->nodes)
- nod->xiExt = nod->xi;
- ddStar[1] = this->runViscous();
- dCdf[1] = vsol->Cdf;
-
- reg->nodes[j]->pos[1] = saveY;
- section->computeLocalCoordinates(reg);
- for (auto &nod : reg->nodes)
- nod->xiExt = nod->xi;
-
- for (int k = 0; k < ddStar[0].size(); ++k)
- T_dStar_x.push_back(Eigen::Triplet<double>(
- k, rowj * isol->pbl->nDim + 1,
- -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
- dCdf_x(rowj * isol->pbl->nDim + 1) += (dCdf[0] - dCdf[1]) / (2. * delta);
- }
- }
- }
- }
- dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end());
- dRdStar_x.prune(0.);
-
- //**************************************************************************//
- // dRblw_x //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- std::deque<Eigen::Triplet<double>> T_blw_x;
- int offSetElms = 0;
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- if (reg->getName() == "wake")
- iReg = 1;
- else if (reg->getName() == "upper" || reg->getName() == "lower")
- iReg = 0;
- else
- throw std::runtime_error("Unknown section name");
- std::vector<std::vector<double>> gradX = reg->evalGradwrtX();
- std::vector<std::vector<double>> gradY = reg->evalGradwrtY();
- for (size_t i = 0; i < gradX.size(); ++i)
- for (size_t j = 0; j < gradX[i].size(); ++j)
- {
- int rowj = isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row;
- T_blw_x.push_back(Eigen::Triplet<double>(
- i + offSetElms, rowj * isol->pbl->nDim + 0, -gradX[i][j]));
- T_blw_x.push_back(Eigen::Triplet<double>(
- i + offSetElms, rowj * isol->pbl->nDim + 1, -gradY[i][j]));
- }
- offSetElms += reg->getnElms();
- }
- }
- }
- dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end());
- dRblw_x.prune(0.);
-
- //**************************************************************************//
- // dRM_x, dRrho_x, dRv_x //
- //--------------------------------------------------------------------------//
- // Analytical //
- //**************************************************************************//
- auto pbl = isol->pbl;
-
- std::deque<Eigen::Triplet<double>> T_dM_x;
- std::deque<Eigen::Triplet<double>> T_drho_x;
- std::deque<Eigen::Triplet<double>> T_dv_x;
-
- int jj = 0;
- iReg = 0;
-
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (auto ® : section->regions)
- {
- if (reg->getName() == "wake")
- iReg = 1;
- else if (reg->getName() == "upper" || reg->getName() == "lower")
- iReg = 0;
- else
- throw std::runtime_error(
- "gradientswrtInviscidFlow: Unknown section name");
- for (size_t iNode = 0; iNode < reg->getnNodes(); ++iNode)
- {
- tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[reg->nodes[iNode]->row];
- double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size());
- for (auto e : pbl->fluid->neMap[srfNode])
- {
- // dv/dx = dv/dgradPhi * dgradPhi/dx = dgradPhi/dx = (-J^(-1) *
- // dJ/dx_j)*grad_xPhi dM/dx = dM/dgradPhi * dgradPhi/dx = dM/dgradPhi *
- // gradPhi * dv/dx drho/dx = drho/dgradPhi * dgradPhi/dx = drho/dgradPhi
- // * gradPhi * dv/dx
- for (size_t inodOnElm = 0; inodOnElm < e->nodes.size(); inodOnElm++)
- for (int idim = 0; idim < nDim; ++idim)
- {
- Eigen::VectorXd gradV = (-e->getJinv(0)) *
- e->getGradJ(0)[inodOnElm * nDim + idim] *
- e->computeGradient(isol->phi, 0);
- double gradM = pbl->fluid->mach->evalGrad(*e, isol->phi, 0) *
- e->computeGradient(isol->phi, 0).transpose() * gradV;
- double gradrho = pbl->fluid->rho->evalGrad(*e, isol->phi, 0) *
- e->computeGradient(isol->phi, 0).transpose() *
- gradV;
- T_dM_x.push_back(Eigen::Triplet<double>(
- jj, e->nodes[inodOnElm]->row * nDim + idim, -gradM / nAdjElms));
- T_drho_x.push_back(Eigen::Triplet<double>(
- jj, e->nodes[inodOnElm]->row * nDim + idim,
- -gradrho / nAdjElms));
- for (int idim2 = 0; idim2 < nDim; idim2++)
- T_dv_x.push_back(Eigen::Triplet<double>(
- jj * nDim + idim2, e->nodes[inodOnElm]->row * nDim + idim,
- -gradV(idim2) / nAdjElms));
- }
- }
- ++jj;
- }
- }
- }
- }
- dRM_x.setFromTriplets(T_dM_x.begin(), T_dM_x.end());
- dRrho_x.setFromTriplets(T_drho_x.begin(), T_drho_x.end());
- dRv_x.setFromTriplets(T_dv_x.begin(), T_dv_x.end());
- dRM_x.prune(0.);
- dRrho_x.prune(0.);
- dRv_x.prune(0.);
+ // int nDim = isol->pbl->nDim;
+
+ // //**************************************************************************//
+ // // dCl_x, dCd_x //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // std::vector<std::vector<double>> dCx_x =
+ // iadjoint->computeGradientCoefficientsMesh();
+ // dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size());
+ // dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
+
+ // //**************************************************************************//
+ // // dRphi_x //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // dRphi_x = iadjoint->getRu_X();
+ // //**************************************************************************//
+ // // dRx_x //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // dRx_x = iadjoint->getRx_X();
+
+ // //**************************************************************************//
+ // // dRdStar_x, dCdf_dx //
+ // //--------------------------------------------------------------------------//
+ // // Central finite-difference //
+ // //--------------------------------------------------------------------------//
+ // // We don't do dRdStar_x = dRdStar_loc * dloc_dx because Cdf = f(xi(x), x) //
+ // // and we would have to compute dCdf/dx and dCdf/dloc. Same for deltaStar. //
+ // // It is more computationally expansive to compute //
+ // // ddeltaStar/dx, ddeltaStar/dloc than ddeltaStar_dx, ddeltaStar_dy //
+ // //**************************************************************************//
+ // std::deque<Eigen::Triplet<double>> T_dStar_x;
+ // double delta = 1e-8;
+ // int iReg = 0;
+ // std::vector<Eigen::VectorXd> ddStar(2,
+ // Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ // std::vector<double> dCdf(2, 0.);
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // if (reg->getName() == "wake")
+ // iReg = 1;
+ // else if (reg->getName() == "upper" || reg->getName() == "lower")
+ // iReg = 0;
+ // else
+ // throw std::runtime_error("Unknown section name");
+ // for (size_t j = 0; j < reg->getnNodes(); ++j)
+ // {
+ // int rowj = isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row;
+ // // x
+ // double saveX = reg->nodes[j]->pos[0];
+
+ // reg->nodes[j]->pos[0] = saveX + delta;
+ // section->computeLocalCoordinates(reg);
+ // for (auto &nod : reg->nodes)
+ // nod->xiExt = nod->xi;
+ // ddStar[0] = this->runViscous();
+ // dCdf[0] = vsol->Cdf;
+ // reg->nodes[j]->pos[0] = saveX - delta;
+ // section->computeLocalCoordinates(reg);
+ // for (auto &nod : reg->nodes)
+ // nod->xiExt = nod->xi;
+ // ddStar[1] = this->runViscous();
+ // dCdf[1] = vsol->Cdf;
+
+ // reg->nodes[j]->pos[0] = saveX;
+ // section->computeLocalCoordinates(reg);
+ // for (auto &nod : reg->nodes)
+ // nod->xiExt = nod->xi;
+
+ // for (int k = 0; k < ddStar[0].size(); ++k)
+ // T_dStar_x.push_back(Eigen::Triplet<double>(
+ // k, rowj * isol->pbl->nDim + 0,
+ // -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
+ // dCdf_x(rowj * isol->pbl->nDim + 0) += (dCdf[0] - dCdf[1]) / (2. * delta);
+ // // y
+ // double saveY = reg->nodes[j]->pos[1];
+
+ // reg->nodes[j]->pos[1] = saveY + delta;
+ // section->computeLocalCoordinates(reg);
+ // for (auto &nod : reg->nodes)
+ // nod->xiExt = nod->xi;
+ // ddStar[0] = this->runViscous();
+ // dCdf[0] = vsol->Cdf;
+ // reg->nodes[j]->pos[1] = saveY - delta;
+ // section->computeLocalCoordinates(reg);
+ // for (auto &nod : reg->nodes)
+ // nod->xiExt = nod->xi;
+ // ddStar[1] = this->runViscous();
+ // dCdf[1] = vsol->Cdf;
+
+ // reg->nodes[j]->pos[1] = saveY;
+ // section->computeLocalCoordinates(reg);
+ // for (auto &nod : reg->nodes)
+ // nod->xiExt = nod->xi;
+
+ // for (int k = 0; k < ddStar[0].size(); ++k)
+ // T_dStar_x.push_back(Eigen::Triplet<double>(
+ // k, rowj * isol->pbl->nDim + 1,
+ // -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
+ // dCdf_x(rowj * isol->pbl->nDim + 1) += (dCdf[0] - dCdf[1]) / (2. * delta);
+ // }
+ // }
+ // }
+ // }
+ // dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end());
+ // dRdStar_x.prune(0.);
+
+ // //**************************************************************************//
+ // // dRblw_x //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // std::deque<Eigen::Triplet<double>> T_blw_x;
+ // int offSetElms = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // if (reg->getName() == "wake")
+ // iReg = 1;
+ // else if (reg->getName() == "upper" || reg->getName() == "lower")
+ // iReg = 0;
+ // else
+ // throw std::runtime_error("Unknown section name");
+ // std::vector<std::vector<double>> gradX = reg->evalGradwrtX();
+ // std::vector<std::vector<double>> gradY = reg->evalGradwrtY();
+ // for (size_t i = 0; i < gradX.size(); ++i)
+ // for (size_t j = 0; j < gradX[i].size(); ++j)
+ // {
+ // int rowj = isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row;
+ // T_blw_x.push_back(Eigen::Triplet<double>(
+ // i + offSetElms, rowj * isol->pbl->nDim + 0, -gradX[i][j]));
+ // T_blw_x.push_back(Eigen::Triplet<double>(
+ // i + offSetElms, rowj * isol->pbl->nDim + 1, -gradY[i][j]));
+ // }
+ // offSetElms += reg->getnElms();
+ // }
+ // }
+ // }
+ // dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end());
+ // dRblw_x.prune(0.);
+
+ // //**************************************************************************//
+ // // dRM_x, dRrho_x, dRv_x //
+ // //--------------------------------------------------------------------------//
+ // // Analytical //
+ // //**************************************************************************//
+ // auto pbl = isol->pbl;
+
+ // std::deque<Eigen::Triplet<double>> T_dM_x;
+ // std::deque<Eigen::Triplet<double>> T_drho_x;
+ // std::deque<Eigen::Triplet<double>> T_dv_x;
+
+ // int jj = 0;
+ // iReg = 0;
+
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ® : section->regions)
+ // {
+ // if (reg->getName() == "wake")
+ // iReg = 1;
+ // else if (reg->getName() == "upper" || reg->getName() == "lower")
+ // iReg = 0;
+ // else
+ // throw std::runtime_error(
+ // "gradientswrtInviscidFlow: Unknown section name");
+ // for (size_t iNode = 0; iNode < reg->getnNodes(); ++iNode)
+ // {
+ // tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[reg->nodes[iNode]->row];
+ // double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ // for (auto e : pbl->fluid->neMap[srfNode])
+ // {
+ // // dv/dx = dv/dgradPhi * dgradPhi/dx = dgradPhi/dx = (-J^(-1) *
+ // // dJ/dx_j)*grad_xPhi dM/dx = dM/dgradPhi * dgradPhi/dx = dM/dgradPhi *
+ // // gradPhi * dv/dx drho/dx = drho/dgradPhi * dgradPhi/dx = drho/dgradPhi
+ // // * gradPhi * dv/dx
+ // for (size_t inodOnElm = 0; inodOnElm < e->nodes.size(); inodOnElm++)
+ // for (int idim = 0; idim < nDim; ++idim)
+ // {
+ // Eigen::VectorXd gradV = (-e->getJinv(0)) *
+ // e->getGradJ(0)[inodOnElm * nDim + idim] *
+ // e->computeGradient(isol->phi, 0);
+ // double gradM = pbl->fluid->mach->evalGrad(*e, isol->phi, 0) *
+ // e->computeGradient(isol->phi, 0).transpose() * gradV;
+ // double gradrho = pbl->fluid->rho->evalGrad(*e, isol->phi, 0) *
+ // e->computeGradient(isol->phi, 0).transpose() *
+ // gradV;
+ // T_dM_x.push_back(Eigen::Triplet<double>(
+ // jj, e->nodes[inodOnElm]->row * nDim + idim, -gradM / nAdjElms));
+ // T_drho_x.push_back(Eigen::Triplet<double>(
+ // jj, e->nodes[inodOnElm]->row * nDim + idim,
+ // -gradrho / nAdjElms));
+ // for (int idim2 = 0; idim2 < nDim; idim2++)
+ // T_dv_x.push_back(Eigen::Triplet<double>(
+ // jj * nDim + idim2, e->nodes[inodOnElm]->row * nDim + idim,
+ // -gradV(idim2) / nAdjElms));
+ // }
+ // }
+ // ++jj;
+ // }
+ // }
+ // }
+ // }
+ // dRM_x.setFromTriplets(T_dM_x.begin(), T_dM_x.end());
+ // dRrho_x.setFromTriplets(T_drho_x.begin(), T_drho_x.end());
+ // dRv_x.setFromTriplets(T_dv_x.begin(), T_dv_x.end());
+ // dRM_x.prune(0.);
+ // dRrho_x.prune(0.);
+ // dRv_x.prune(0.);
}
/**
@@ -1241,41 +1241,41 @@ void CoupledAdjoint::gradientwrtMesh()
*/
Eigen::VectorXd CoupledAdjoint::runViscous()
{
- int exitCode = vsol->run();
- if (exitCode != 0 && vsol->verbose > 0)
- std::cout << "vsol terminated with exit code " << exitCode << std::endl;
-
- // Extract deltaStar
- std::vector<Eigen::VectorXd> dStar_p;
-
- // Extract parts in a loop
- for (auto &body : vsol->bodies)
- {
- for (auto §ion : body->sections)
- {
- for (size_t i = 0; i < section->regions.size(); ++i)
- {
- std::vector<double> deltaStar_vec = section->regions[i]->deltaStar;
- Eigen::VectorXd deltaStar_eigen =
- Eigen::Map<Eigen::VectorXd>(deltaStar_vec.data(), deltaStar_vec.size());
- dStar_p.push_back(deltaStar_eigen);
- }
- }
- }
- // Calculate the total size
- int nNs = 0;
- for (const auto &p : dStar_p)
- nNs += p.size();
-
- // Concatenate the vectors
- Eigen::VectorXd deltaStar(nNs);
- int idx = 0;
- for (const auto &p : dStar_p)
- {
- deltaStar.segment(idx, p.size()) = p;
- idx += p.size();
- }
- return deltaStar;
+ // int exitCode = vsol->run();
+ // if (exitCode != 0 && vsol->verbose > 0)
+ // std::cout << "vsol terminated with exit code " << exitCode << std::endl;
+
+ // // Extract deltaStar
+ // std::vector<Eigen::VectorXd> dStar_p;
+
+ // // Extract parts in a loop
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (size_t i = 0; i < section->regions.size(); ++i)
+ // {
+ // std::vector<double> deltaStar_vec = section->regions[i]->deltaStar;
+ // Eigen::VectorXd deltaStar_eigen =
+ // Eigen::Map<Eigen::VectorXd>(deltaStar_vec.data(), deltaStar_vec.size());
+ // dStar_p.push_back(deltaStar_eigen);
+ // }
+ // }
+ // }
+ // // Calculate the total size
+ // int nNs = 0;
+ // for (const auto &p : dStar_p)
+ // nNs += p.size();
+
+ // // Concatenate the vectors
+ // Eigen::VectorXd deltaStar(nNs);
+ // int idx = 0;
+ // for (const auto &p : dStar_p)
+ // {
+ // deltaStar.segment(idx, p.size()) = p;
+ // idx += p.size();
+ // }
+ // return deltaStar;
}
/**
diff --git a/blast/src/blDriver.cpp b/blast/src/blDriver.cpp
index a23f3dc416bb9a5f71b1e40335156a9c9fcde417..ba33f166d20497813e3c07f486c72ec65fff3c29 100644
--- a/blast/src/blDriver.cpp
+++ b/blast/src/blDriver.cpp
@@ -1,3 +1,4 @@
+#include "blast.h"
#include "blDriver.h"
#include "blBoundaryLayer.h"
#include "blSolver.h"
@@ -64,7 +65,7 @@ int Driver::run()
// unsuccessful, failed; >0: unsuccessful nb iter).
int numberFailedPerBody = 0;
- double maxMach = 0.;
+ bldouble maxMach = 0.;
solverExitCode = 0;
@@ -109,7 +110,7 @@ int Driver::run()
for (auto &nod : reg->nodes)
nod->setRegime(1);
// Impose wake boundary condition.
- std::vector<std::vector<double>> teConditions;
+ std::vector<std::vector<bldouble>> teConditions;
for (auto &iregion : section->regions)
if (iregion->getName() == "upper" || iregion->getName() == "lower")
teConditions.push_back(iregion->getTeConditions());
@@ -200,7 +201,7 @@ int Driver::run()
void Driver::checkWaves(size_t inod, std::shared_ptr<BoundaryLayer> &bl)
{
- double logRet_crit =
+ bldouble logRet_crit =
2.492 * std::pow(1. / (bl->Hk[inod] - 1), 0.43) +
0.7 * (std::tanh(14 * (1. / (bl->Hk[inod] - 1)) - 9.24) + 1.0);
@@ -229,7 +230,7 @@ void Driver::averageTransition(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
size_t nVar = bl->getnVar();
// Save laminar solution.
- std::vector<double> lamSol(nVar, 0.0);
+ std::vector<bldouble> lamSol(nVar, 0.0);
for (size_t k = 0; k < lamSol.size(); ++k)
lamSol[k] = bl->u[inod * nVar + k];
@@ -245,16 +246,16 @@ void Driver::averageTransition(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
(bl->u[inod * nVar + 2] - bl->u[(inod - 1) * nVar + 2])) +
bl->nodes[inod - 1]->xoc);
- double avTurb =
+ bldouble avTurb =
(bl->nodes[inod]->xoc - bl->getxtr()) / (bl->nodes[inod]->xoc - bl->nodes[inod - 1]->xoc);
if (avTurb < 0. - 1e-12)
avTurb = 0.;
if (avTurb > 1. + 1e-12)
avTurb = 1.;
- double avLam = 1. - avTurb;
+ bldouble avLam = 1. - avTurb;
// Impose boundary condition.
- double Cteq_trans;
+ bldouble Cteq_trans;
bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(inod - 1) * nVar + 0],
bl->u[(inod - 1) * nVar + 1], 0.03,
bl->Ret[inod - 1], bl->Me[inod - 1],
@@ -295,7 +296,7 @@ void Driver::averageTransition(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
// Recompute closures. (The turbulent BC @ inod - 1 does not influence
// laminar closure relations).
- std::vector<double> lamParam(8, 0.);
+ std::vector<bldouble> lamParam(8, 0.);
bl->closSolver->laminarClosures(
lamParam, bl->u[(inod - 1) * nVar + 0], bl->u[(inod - 1) * nVar + 1],
bl->Ret[inod - 1], bl->Me[inod - 1], bl->getName());
@@ -308,7 +309,7 @@ void Driver::averageTransition(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
bl->ctEq[inod - 1] = lamParam[6];
bl->us[inod - 1] = lamParam[7];
- std::vector<double> turbParam(8, 0.);
+ std::vector<bldouble> turbParam(8, 0.);
bl->closSolver->turbulentClosures(
turbParam, bl->u[(inod)*nVar + 0], bl->u[(inod)*nVar + 1],
bl->u[(inod)*nVar + 4], bl->Ret[inod], bl->Me[inod], bl->getName());
diff --git a/blast/src/blDriver.h b/blast/src/blDriver.h
index 7c5ffd70ae2d99eb2d36eeaaf856316e09473c42..79356e9e8ba26688c3046522ed66ba6d6c00e2eb 100644
--- a/blast/src/blDriver.h
+++ b/blast/src/blDriver.h
@@ -23,18 +23,16 @@ public:
std::vector<std::shared_ptr<Body>> bodies; ///< Boundary layer regions sorted by body and section.
unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1.
- double Cdt;
- double Cdf;
- double Cdp;
-
private:
- double Re; ///< Freestream Reynolds number.
- double CFL0; ///< Initial CFL number for pseudo time integration.
- std::vector<double> spans; ///< Wing Span (Used only for drag computation, not used if 2D
- ///< case).
+ bldouble Re; ///< Freestream Reynolds number.
+ bldouble CFL0; ///< Initial CFL number for pseudo time integration.
Solver *tSolver; ///< Pseudo-time solver.
int solverExitCode; ///< Exit code of viscous calculations.
+ bldouble Cdt;
+ bldouble Cdf;
+ bldouble Cdp;
+
protected:
fwk::Timers tms; ///< Internal timers
@@ -47,8 +45,11 @@ public:
void add(std::shared_ptr<Body> body);
// Getters.
- double getRe() const { return Re; }
- double getMinf() const { return tSolver->getMinf(); }
+ bldouble getRe() const { return Re; }
+ bldouble getMinf() const { return tSolver->getMinf(); }
+ bldouble getCdf() const { return Cdf; }
+ bldouble getCdt() const { return Cdt; }
+ bldouble getCdp() const { return Cdp; }
// Setters.
void setVerbose(unsigned int _verbose) { verbose = _verbose; }
diff --git a/blast/src/blFluxes.cpp b/blast/src/blFluxes.cpp
index 2578f52d04862edd5428d6bef9692ac051b1ac56..8b974ba90330a5bebac1b8a533a4b18c2bcef64d 100644
--- a/blast/src/blFluxes.cpp
+++ b/blast/src/blFluxes.cpp
@@ -1,53 +1,87 @@
+#include "blast.h"
#include "blFluxes.h"
#include "blBoundaryLayer.h"
#include "blNode.h"
-#include <Eigen/Dense>
+
+#ifdef USE_CODI
+#include "codi.hpp"
+#endif
using namespace blast;
-using namespace Eigen;
/**
* @brief Compute residual vector of the integral boundary layer equations.
*
* @param inod Marker id.
* @param bl Boundary layer region.
- * @return VectorXd
+ * @return blVectorXd
*/
-VectorXd Fluxes::computeResiduals(size_t inod, std::shared_ptr<BoundaryLayer> &bl)
+blVectorXd Fluxes::computeResiduals(size_t inod, std::shared_ptr<BoundaryLayer> &bl)
{
size_t nVar = bl->getnVar();
- std::vector<double> u(nVar, 0.);
+ std::vector<bldouble> u(nVar, 0.);
for (size_t k = 0; k < nVar; ++k)
u[k] = bl->u[inod * nVar + k];
return blLaws(inod, bl, u);
}
/**
- * @brief Compute Jacobian of the integral boundary layer equations.
+ * @brief Compute Jacobian of the integral boundary layer equations using
+ * automatic differentiation or finite differences.
*
* @param inod Marker id.
* @param bl Boundary layer region.
* @param sysRes Residual of the IBL at point inod.
* @param eps Pertubation from current solution used to compute the Jacobian.
- * @return MatrixXd
+ * @return blMatrixXd
*/
-MatrixXd Fluxes::computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
- VectorXd const &sysRes, double eps)
+blMatrixXd Fluxes::computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
+ blVectorXd const &sysRes, bldouble eps)
{
size_t nVar = bl->getnVar();
- MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
- std::vector<double> uUp(nVar, 0.);
+ blMatrixXd JacMatrix = blMatrixXd::Zero(nVar, nVar);
+ std::vector<bldouble> sol(nVar, 0.);
for (size_t k = 0; k < nVar; ++k)
- uUp[k] = bl->u[inod * nVar + k];
+ sol[k] = bl->u[inod * nVar + k];
- double varSave = 0.;
- for (size_t iVar = 0; iVar < nVar; ++iVar)
- {
- varSave = uUp[iVar];
- uUp[iVar] += eps;
- JacMatrix.col(iVar) = (blLaws(inod, bl, uUp) - sysRes) / eps;
- uUp[iVar] = varSave;
- }
+ #ifdef USE_CODI
+ // Compute the Jacobian using automatic differentiation.
+ bltape &tape = bldouble::getTape();
+ tape.reset();
+ tape.setActive();
+ // Register inputs
+ for (auto &ui : sol)
+ tape.registerInput(ui);
+ // Compute residuals
+ blVectorXd res = blLaws(inod, bl, sol);
+ // Register outputs
+ for (auto& ri : res)
+ tape.registerOutput(ri);
+ // Desactivate tape
+ tape.setPassive();
+
+ // Form the Jacobian
+ for (size_t i = 0; i < 5; ++i)
+ {
+ res[i].gradient() = 1.0;
+
+ tape.evaluate();
+ for(size_t j = 0; j < 5; ++j) {
+ JacMatrix(i, j) = sol[j].getGradient();
+ }
+ tape.clearAdjoints();
+ }
+ #else
+ // Compute the Jacobian using finite differences.
+ bldouble varSave = 0.;
+ for (size_t iVar = 0; iVar < nVar; ++iVar)
+ {
+ varSave = sol[iVar];
+ sol[iVar] += eps;
+ JacMatrix.col(iVar) = (blLaws(inod, bl, sol) - sysRes) / eps;
+ sol[iVar] = varSave;
+ }
+ #endif
return JacMatrix;
}
@@ -57,17 +91,17 @@ MatrixXd Fluxes::computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &bl
* @param inod Node number.
* @param bl Boundary layer region.
* @param u Solution vector at point inod.
- * @return VectorXd
+ * @return blVectorXd
*/
-VectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
- std::vector<double> u)
+blVectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
+ std::vector<bldouble> &u)
{
size_t nVar = bl->getnVar();
auto &currNode = bl->nodes[inod];
auto &prevNode = bl->nodes[inod - 1];
// Dissipation ratio.
- double dissipRatio = 0.;
+ bldouble dissipRatio = 0.;
if (bl->getName() == "wake")
dissipRatio = 0.9;
else
@@ -81,7 +115,7 @@ VectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
if (currNode->getRegime() == 0)
{
// Laminar closure relations.
- std::vector<double> lamParam(8, 0.);
+ std::vector<bldouble> lamParam(8, 0.);
bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[inod],
bl->Me[inod], bl->getName());
bl->Hstar[inod] = lamParam[0];
@@ -96,7 +130,7 @@ VectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
else if (currNode->getRegime() == 1)
{
// Turbulent closure relations.
- std::vector<double> turbParam(8, 0.);
+ std::vector<bldouble> turbParam(8, 0.);
bl->closSolver->turbulentClosures(
turbParam, u[0], u[1], u[4], bl->Ret[inod], bl->Me[inod], bl->getName());
bl->Hstar[inod] = turbParam[0];
@@ -112,24 +146,24 @@ VectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
throw std::runtime_error("Wrong regime");
// Gradients computation.
- double dx = (currNode->xi - prevNode->xi);
- double dxExt = (currNode->xiExt - prevNode->xiExt);
- double dt_dx = (u[0] - bl->u[(inod - 1) * nVar + 0]) / dx;
- double dH_dx = (u[1] - bl->u[(inod - 1) * nVar + 1]) / dx;
- double due_dx = (u[3] - bl->u[(inod - 1) * nVar + 3]) / dx;
- double dct_dx = (u[4] - bl->u[(inod - 1) * nVar + 4]) / dx;
- double dHstar_dx = (bl->Hstar[inod] - bl->Hstar[inod - 1]) / dx;
+ bldouble dx = (currNode->xi - prevNode->xi);
+ bldouble dxExt = (currNode->xiExt - prevNode->xiExt);
+ bldouble dt_dx = (u[0] - bl->u[(inod - 1) * nVar + 0]) / dx;
+ bldouble dH_dx = (u[1] - bl->u[(inod - 1) * nVar + 1]) / dx;
+ bldouble due_dx = (u[3] - bl->u[(inod - 1) * nVar + 3]) / dx;
+ bldouble dct_dx = (u[4] - bl->u[(inod - 1) * nVar + 4]) / dx;
+ bldouble dHstar_dx = (bl->Hstar[inod] - bl->Hstar[inod - 1]) / dx;
- double dueExt_dx = (bl->vt[inod] - bl->vt[inod - 1]) / dxExt;
- double ddeltaStarExt_dx =
+ bldouble dueExt_dx = (bl->vt[inod] - bl->vt[inod - 1]) / dxExt;
+ bldouble ddeltaStarExt_dx =
(bl->deltaStarExt[inod] - bl->deltaStarExt[inod - 1]) / dxExt;
- double c = 4. / (M_PI * dx);
- double cExt = 4. / (M_PI * dxExt);
+ bldouble c = 4. / (M_PI * dx);
+ bldouble cExt = 4. / (M_PI * dxExt);
// Amplification routine.
- double dN_dx = 0.0;
- double ax = 0.0;
+ bldouble dN_dx = 0.0;
+ bldouble ax = 0.0;
if (bl->getxtrF() == -1.0 && currNode->getRegime() == 0)
{
// No forced transition and laminar regime.
@@ -137,15 +171,15 @@ VectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
dN_dx = (u[2] - bl->u[(inod - 1) * nVar + 2]) / dx;
}
- double Mea = bl->Me[inod];
- double Hstara = bl->Hstar[inod];
- double Hstar2a = bl->Hstar2[inod];
- double Hka = bl->Hk[inod];
- double cfa = bl->cf[inod];
- double cda = bl->cd[inod];
- double deltaa = bl->delta[inod];
- double ctEqa = bl->ctEq[inod];
- double usa = bl->us[inod];
+ bldouble Mea = bl->Me[inod];
+ bldouble Hstara = bl->Hstar[inod];
+ bldouble Hstar2a = bl->Hstar2[inod];
+ bldouble Hka = bl->Hk[inod];
+ bldouble cfa = bl->cf[inod];
+ bldouble cda = bl->cd[inod];
+ bldouble deltaa = bl->delta[inod];
+ bldouble ctEqa = bl->ctEq[inod];
+ bldouble usa = bl->us[inod];
//--------------------------------------------------------------------------------//
// Integral boundary layer equations. //
@@ -187,7 +221,7 @@ VectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
// timeMatrix(4, 4) = 1.0;
//--------------------------------------------------------------------------------//
- VectorXd result = VectorXd::Zero(5);
+ blVectorXd result = blVectorXd::Zero(5);
if (currNode->getRegime() == 0)
{
@@ -298,28 +332,28 @@ VectorXd Fluxes::blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
* @param Hk Kinematic shape parameter.
* @param Ret Reynolds number based on the momentum thickness.
* @param theta Momentum thickness.
- * @return double
+ * @return bldouble
*/
-double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const
+bldouble Fluxes::amplificationRoutine(bldouble Hk, bldouble Ret, bldouble theta) const
{
- double dgr = 0.08;
- double Hk1 = 3.5;
- double Hk2 = 4.;
- double Hmi = (1 / (Hk - 1));
- double logRet_crit =
+ bldouble dgr = 0.08;
+ bldouble Hk1 = 3.5;
+ bldouble Hk2 = 4.;
+ bldouble Hmi = (1 / (Hk - 1));
+ bldouble logRet_crit =
2.492 * std::pow(1. / (Hk - 1.), 0.43) +
0.7 * (std::tanh(14. * Hmi - 9.24) +
1.0); // value at which the amplification starts to grow
- double ax = 0.;
+ bldouble ax = 0.;
if (Ret <= 0.)
Ret = 1.;
if (std::log10(Ret) < logRet_crit - dgr)
ax = 0.;
else
{
- double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2. * dgr);
- double Rfac = 0.;
+ bldouble rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2. * dgr);
+ bldouble Rfac = 0.;
if (rNorm <= 0)
Rfac = 0.;
if (rNorm >= 1)
@@ -328,25 +362,25 @@ double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const
Rfac = 3. * rNorm * rNorm - 2. * rNorm * rNorm * rNorm;
// Normal envelope amplification rate
- double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi +
+ bldouble m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi +
0.1 * std::exp(-20. * Hmi);
- double arg = 3.87 * Hmi - 2.52;
- double dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg));
+ bldouble arg = 3.87 * Hmi - 2.52;
+ bldouble dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg));
ax = (m * dndRet / theta) * Rfac;
// Set correction for separated profiles
if (Hk > Hk1)
{
- double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
- double Hfac = 0.;
+ bldouble Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
+ bldouble Hfac = 0.;
if (Hnorm >= 1.)
Hfac = 1.;
else
Hfac = 3. * Hnorm * Hnorm - 2. * Hnorm * Hnorm * Hnorm;
- double ax1 = ax;
- double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5.));
- double tnr = std::tanh(1.2 * (std::log10(Ret) - gro));
- double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1.), 1.5)) / theta;
+ bldouble ax1 = ax;
+ bldouble gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5.));
+ bldouble tnr = std::tanh(1.2 * (std::log10(Ret) - gro));
+ bldouble ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1.), 1.5)) / theta;
if (ax2 < 0.)
ax2 = 0.;
ax = Hfac * ax2 + (1. - Hfac) * ax1;
diff --git a/blast/src/blFluxes.h b/blast/src/blFluxes.h
index 4e8d5b12932ed1d276fefa9975ab5ee3cdff4e1c..dee71ab5b46e1ad4101ed88136bf2cb1c33ca83b 100644
--- a/blast/src/blFluxes.h
+++ b/blast/src/blFluxes.h
@@ -3,8 +3,6 @@
#include "blBoundaryLayer.h"
#include "blast.h"
-#include <Eigen/Dense>
-
#include <memory>
namespace blast
@@ -18,19 +16,19 @@ namespace blast
class BLAST_API Fluxes
{
public:
- double Re; ///< Freestream Reynolds number.
+ bldouble Re; ///< Freestream Reynolds number.
public:
- Fluxes(double _Re) : Re(_Re) {}
+ Fluxes(bldouble _Re) : Re(_Re) {}
~Fluxes() {}
- Eigen::VectorXd computeResiduals(size_t inod, std::shared_ptr<BoundaryLayer> &bl);
- Eigen::MatrixXd computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
- Eigen::VectorXd const &sysRes, double eps);
+ blVectorXd computeResiduals(size_t inod, std::shared_ptr<BoundaryLayer> &bl);
+ blMatrixXd computeJacobian(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
+ blVectorXd const &sysRes, bldouble eps);
private:
- Eigen::VectorXd blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
- std::vector<double> u);
- double amplificationRoutine(double Hk, double Ret, double theta) const;
+ blVectorXd blLaws(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
+ std::vector<bldouble> &u);
+ bldouble amplificationRoutine(bldouble Hk, bldouble Ret, bldouble theta) const;
};
} // namespace blast
#endif // BLFLUXES_H
diff --git a/blast/src/blNode.cpp b/blast/src/blNode.cpp
index dd110d525a84394c91fd5761a2111f7bc58ea722..2628aa78330694dd33999e7b51f86cca19ce0771 100644
--- a/blast/src/blNode.cpp
+++ b/blast/src/blNode.cpp
@@ -16,6 +16,7 @@
// Adapted from tbox wNode by Romain Boman.
+#include "blast.h"
#include "blNode.h"
#include <Eigen/Dense>
@@ -26,7 +27,7 @@ void blNode::write(std::ostream &out) const
out << "node #" << no << " = " << pos.transpose() << " row=" << row;
}
-blNode::blNode(int n, Eigen::Vector3d const &p, int r) : wObject(), no(n), pos(p), row(r)
+blNode::blNode(int n, blVector3d const &p, int r) : wObject(), no(n), pos(p), row(r)
{
xoc = 0.0;
xi = 0.0;
diff --git a/blast/src/blNode.h b/blast/src/blNode.h
index 026c068acf6f24ccd0b4b305132baf8c4d5f997c..84d16040895ce450202b9fda3b483d739a0f4884 100644
--- a/blast/src/blNode.h
+++ b/blast/src/blNode.h
@@ -24,7 +24,6 @@
#include <vector>
#include <algorithm>
#include <memory>
-#include <Eigen/Dense>
namespace blast
{
@@ -36,18 +35,18 @@ class BLAST_API blNode : public fwk::wObject
{
public:
int no;
- Eigen::Vector3d pos;
+ blVector3d pos;
int row;
- double xoc; ///< x/c coordinate in the global frame of reference.
- double xi; ///< coordinate in the local frame of reference.
- double xiExt; ///< coordinate in the local frame of reference, at the edge of the boundary layer.
+ bldouble xoc; ///< x/c coordinate in the global frame of reference.
+ bldouble xi; ///< coordinate in the local frame of reference.
+ bldouble xiExt; ///< coordinate in the local frame of reference, at the edge of the boundary layer.
private:
int regime; ///< Regime of the boundary layer (0=laminar, 1=turbulent).
public:
- blNode(int n, Eigen::Vector3d const &p, int r);
+ blNode(int n, blVector3d const &p, int r);
void setRegime(int const r) { regime = r; }
int getRegime() const { return regime; }
diff --git a/blast/src/blSection.cpp b/blast/src/blSection.cpp
index bcae644e5a8495f764c458432e646ae918eca142..3de8e6dc422ab18b5609e1848427808b33b9f237 100644
--- a/blast/src/blSection.cpp
+++ b/blast/src/blSection.cpp
@@ -14,9 +14,9 @@
* limitations under the License.
*/
+#include "blast.h"
#include "blSection.h"
#include "blBoundaryLayer.h"
-#include <Eigen/Dense>
using namespace blast;
@@ -32,8 +32,8 @@ void Section::add(std::shared_ptr<BoundaryLayer> bl)
void Section::computeLocalCoordinates()
{
// Compute minimum x and local chord
- double xmin = 0.0;
- double xmax = 0.0;
+ bldouble xmin = 0.0;
+ bldouble xmax = 0.0;
for (auto &bl : regions)
{
if (bl->getName() == "wake")
@@ -46,7 +46,7 @@ void Section::computeLocalCoordinates()
xmax = node->pos[0];
}
}
- double chord = xmax - xmin;
+ bldouble chord = xmax - xmin;
if (chord <= 0.)
throw std::runtime_error("Chord length is zero or negative.");
@@ -58,7 +58,7 @@ void Section::computeLocalCoordinates()
{
auto &nod1 = bl->nodes[ielm + 1];
auto &nod2 = bl->nodes[ielm];
- double dx = std::sqrt((nod1->pos[0] - nod2->pos[0]) * (nod1->pos[0] - nod2->pos[0]) + (nod1->pos[1] - nod2->pos[1]) * (nod1->pos[1] - nod2->pos[1]));
+ bldouble dx = std::sqrt((nod1->pos[0] - nod2->pos[0]) * (nod1->pos[0] - nod2->pos[0]) + (nod1->pos[1] - nod2->pos[1]) * (nod1->pos[1] - nod2->pos[1]));
nod1->xi = nod2->xi + dx;
nod1->xiExt = nod1->xi;
}
@@ -71,8 +71,8 @@ void Section::computeLocalCoordinates()
void Section::computeLocalCoordinates(std::shared_ptr<BoundaryLayer> ®)
{
// Compute minimum x and local chord
- double xmin = 0.0;
- double xmax = 0.0;
+ bldouble xmin = 0.0;
+ bldouble xmax = 0.0;
for (auto &bl : regions)
{
if (bl->getName() == "wake")
@@ -85,7 +85,7 @@ void Section::computeLocalCoordinates(std::shared_ptr<BoundaryLayer> ®)
xmax = node->pos[0];
}
}
- double chord = xmax - xmin;
+ bldouble chord = xmax - xmin;
if (chord <= 0.)
throw std::runtime_error("Chord length is zero or negative.");
@@ -95,7 +95,7 @@ void Section::computeLocalCoordinates(std::shared_ptr<BoundaryLayer> ®)
{
auto &nod1 = reg->nodes[ielm + 1];
auto &nod2 = reg->nodes[ielm];
- double dx = std::sqrt((nod1->pos[0] - nod2->pos[0]) * (nod1->pos[0] - nod2->pos[0]) + (nod1->pos[1] - nod2->pos[1]) * (nod1->pos[1] - nod2->pos[1]));
+ bldouble dx = std::sqrt((nod1->pos[0] - nod2->pos[0]) * (nod1->pos[0] - nod2->pos[0]) + (nod1->pos[1] - nod2->pos[1]) * (nod1->pos[1] - nod2->pos[1]));
nod1->xi = nod2->xi + dx;
// nod1->xiExt = nod1->xi;
}
@@ -122,7 +122,7 @@ void Section::computeDrag()
Cdf = 0.0;
for (auto ® : regions)
{
- std::vector<double> cf_side(reg->getnNodes(), 0.0);
+ std::vector<bldouble> cf_side(reg->getnNodes(), 0.0);
for (size_t k = 0; k < reg->getnNodes(); ++k)
cf_side[k] = reg->vt[k] * reg->vt[k] * reg->cf[k];
// Integrate cf_side over reg->nodes[i]->pos[0] and pos[1]
diff --git a/blast/src/blSection.h b/blast/src/blSection.h
index 5039f4902ad2544f72abff03b942c3c922ca19c1..0c8c7ea8fd35c0d82478005751dc08404e63eb0c 100644
--- a/blast/src/blSection.h
+++ b/blast/src/blSection.h
@@ -19,9 +19,9 @@ class BLAST_API Section : public fwk::wSharedObject
{
private:
- double Cdt = 0.0; ///< Total drag coefficient.
- double Cdf = 0.0; ///< Friction drag coefficient.
- double Cdp = 0.0; ///< Pressure drag coefficient.
+ bldouble Cdt = 0.0; ///< Total drag coefficient.
+ bldouble Cdf = 0.0; ///< Friction drag coefficient.
+ bldouble Cdp = 0.0; ///< Pressure drag coefficient.
public:
Section();
@@ -32,9 +32,9 @@ public:
void computeLocalCoordinates();
void computeLocalCoordinates(std::shared_ptr<BoundaryLayer> ®);
void computeDrag();
- double getCdt() const { return Cdt; }
- double getCdf() const { return Cdf; }
- double getCdp() const { return Cdp; }
+ bldouble getCdt() const { return Cdt; }
+ bldouble getCdf() const { return Cdf; }
+ bldouble getCdp() const { return Cdp; }
int getRegionIdByName(std::string name) const;
bool hasWake() const;
};
diff --git a/blast/src/blSolver.cpp b/blast/src/blSolver.cpp
index c34e7c02ddc735efa2af561a9dcc42bb40d83568..f0fb70dbad074c56f34db3ab10bd58720ba759e2 100644
--- a/blast/src/blSolver.cpp
+++ b/blast/src/blSolver.cpp
@@ -1,10 +1,9 @@
+#include "blast.h"
#include "blSolver.h"
#include "blBoundaryLayer.h"
#include "blFluxes.h"
-#include <Eigen/Dense>
#include <iostream>
-using namespace Eigen;
using namespace tbox;
using namespace blast;
@@ -19,8 +18,8 @@ using namespace blast;
* @param _itFrozenJac Number of iterations between which the Jacobian is
* frozen.
*/
-Solver::Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt,
- double _tol, unsigned int _itFrozenJac)
+Solver::Solver(bldouble _CFL0, bldouble _Minf, bldouble Re, unsigned int _maxIt,
+ bldouble _tol, unsigned int _itFrozenJac)
{
CFL0 = _CFL0;
maxIt = _maxIt;
@@ -81,27 +80,27 @@ int Solver::integration(size_t inod, std::shared_ptr<BoundaryLayer> &bl)
auto prevNode = bl->nodes[inod - 1];
// Save initial condition.
- std::vector<double> sln0(nVar, 0.);
+ std::vector<bldouble> sln0(nVar, 0.);
for (size_t i = 0; i < nVar; ++i)
sln0[i] = bl->u[inod * nVar + i];
// Initialize time integration variables.
- double dx = currNode->xi - prevNode->xi;
+ bldouble dx = currNode->xi - prevNode->xi;
// Initial time step.
- double CFL = CFL0;
+ bldouble CFL = CFL0;
unsigned int itFrozenJac = itFrozenJac0;
- double timeStep = setTimeStep(CFL, dx, bl->vt[inod]);
- MatrixXd IMatrix(5, 5);
- IMatrix = MatrixXd::Identity(5, 5) / timeStep;
+ bldouble timeStep = setTimeStep(CFL, dx, bl->vt[inod]);
+ blMatrixXd IMatrix(5, 5);
+ IMatrix = blMatrixXd::Identity(5, 5) / timeStep;
// Initial system residuals.
- VectorXd sysRes = sysMatrix->computeResiduals(inod, bl);
- double normSysRes0 = sysRes.norm();
- double normSysRes = normSysRes0;
+ blVectorXd sysRes = sysMatrix->computeResiduals(inod, bl);
+ bldouble normSysRes0 = sysRes.norm();
+ bldouble normSysRes = normSysRes0;
- MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
- VectorXd slnIncr = VectorXd::Zero(5);
+ blMatrixXd JacMatrix = blMatrixXd::Zero(5, 5);
+ blVectorXd slnIncr = blVectorXd::Zero(5);
unsigned int innerIters = 0; // Inner (non-linear) iterations
@@ -129,7 +128,7 @@ int Solver::integration(size_t inod, std::shared_ptr<BoundaryLayer> &bl)
// CFL adaptation.
CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1);
timeStep = setTimeStep(CFL, dx, bl->vt[inod]);
- IMatrix = MatrixXd::Identity(5, 5) / timeStep;
+ IMatrix = blMatrixXd::Identity(5, 5) / timeStep;
innerIters++;
}
@@ -149,16 +148,16 @@ int Solver::integration(size_t inod, std::shared_ptr<BoundaryLayer> &bl)
* @param CFL CFL number.
* @param dx Cell size.
* @param invVel Inviscid velocity.
- * @return double
+ * @return bldouble
*/
-double Solver::setTimeStep(double CFL, double dx, double invVel) const
+bldouble Solver::setTimeStep(bldouble CFL, bldouble dx, bldouble invVel) const
{
// Speed of sound.
- double gradU2 = (invVel * invVel);
- double soundSpeed = computeSoundSpeed(gradU2);
+ bldouble gradU2 = (invVel * invVel);
+ bldouble soundSpeed = computeSoundSpeed(gradU2);
// Velocity of the fastest traveling wave.
- double advVel = soundSpeed + invVel;
+ bldouble advVel = soundSpeed + invVel;
// Time step computation.
return CFL * dx / advVel;
@@ -168,9 +167,9 @@ double Solver::setTimeStep(double CFL, double dx, double invVel) const
* @brief Compute the speed of sound in a cell.
*
* @param gradU2 Inviscid velocity squared in the cell.
- * @return double
+ * @return bldouble
*/
-double Solver::computeSoundSpeed(double gradU2) const
+bldouble Solver::computeSoundSpeed(bldouble gradU2) const
{
return sqrt(1. / (Minf * Minf) + 0.5 * (gamma - 1.) * (1. - gradU2));
}
@@ -185,7 +184,7 @@ double Solver::computeSoundSpeed(double gradU2) const
* @param usePrevPoint if 1, use the solution at previous point instead of sln0.
*/
void Solver::resetSolution(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
- std::vector<double> sln0, bool usePrevPoint)
+ std::vector<bldouble> sln0, bool usePrevPoint)
{
size_t nVar = bl->getnVar();
@@ -206,7 +205,7 @@ void Solver::resetSolution(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
if (bl->nodes[inod]->getRegime() == 0)
{
- std::vector<double> lamParam(8, 0.);
+ std::vector<bldouble> lamParam(8, 0.);
bl->closSolver->laminarClosures(lamParam, bl->u[inod * nVar + 0],
bl->u[inod * nVar + 1], bl->Ret[inod],
bl->Me[inod], bl->getName());
@@ -221,7 +220,7 @@ void Solver::resetSolution(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
}
else if (bl->nodes[inod]->getRegime() == 1)
{
- std::vector<double> turbParam(8, 0.);
+ std::vector<bldouble> turbParam(8, 0.);
bl->closSolver->turbulentClosures(
turbParam, bl->u[inod * nVar + 0], bl->u[inod * nVar + 1],
bl->u[inod * nVar + 4], bl->Ret[inod], bl->Me[inod], bl->getName());
@@ -236,7 +235,7 @@ void Solver::resetSolution(size_t inod, std::shared_ptr<BoundaryLayer> &bl,
}
}
-void Solver::setCFL0(double _CFL0)
+void Solver::setCFL0(bldouble _CFL0)
{
if (_CFL0 > 0)
CFL0 = _CFL0;
diff --git a/blast/src/blSolver.h b/blast/src/blSolver.h
index 5e1753e1ceb38fe5fb06adf835b4a45b0bc10f80..b882fb96e9eec7db17d2dd8630da19c9f8dbf70b 100644
--- a/blast/src/blSolver.h
+++ b/blast/src/blSolver.h
@@ -17,31 +17,31 @@ class BLAST_API Solver
{
private:
- double CFL0; ///< Initial CFL number.
+ bldouble CFL0; ///< Initial CFL number.
unsigned int maxIt; ///< Maximum number of iterations.
- double tol; ///< Convergence tolerance.
+ bldouble tol; ///< Convergence tolerance.
unsigned int itFrozenJac0; ///< Number of iterations between which the
///< Jacobian is frozen.
- double const gamma = 1.4; ///< Air heat capacity ratio.
- double Minf; ///< Freestream Mach number.
+ bldouble const gamma = 1.4; ///< Air heat capacity ratio.
+ bldouble Minf; ///< Freestream Mach number.
Fluxes *sysMatrix;
public:
- Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 100,
- double _tol = 1e-8, unsigned int _itFrozenJac = 1);
+ Solver(bldouble _CFL0, bldouble _Minf, bldouble Re, unsigned int _maxIt = 100,
+ bldouble _tol = 1e-8, unsigned int _itFrozenJac = 1);
~Solver();
void initialCondition(size_t inod, std::shared_ptr<BoundaryLayer> &bl);
int integration(size_t inod, std::shared_ptr<BoundaryLayer> &bl);
// Getters.
- double getCFL0() const { return CFL0; }
- double getMinf() const { return Minf; }
+ bldouble getCFL0() const { return CFL0; }
+ bldouble getMinf() const { return Minf; }
unsigned int getmaxIt() const { return maxIt; }
unsigned int getitFrozenJac() const { return itFrozenJac0; }
// Setters.
- void setCFL0(double _CFL0);
+ void setCFL0(bldouble _CFL0);
void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; };
void setitFrozenJac(unsigned int _itFrozenJac)
{
@@ -49,9 +49,9 @@ public:
}
private:
- double setTimeStep(double CFL, double dx, double invVel) const;
- double computeSoundSpeed(double gradU2) const;
- void resetSolution(size_t inod, std::shared_ptr<BoundaryLayer> &bl, std::vector<double> Sln0,
+ bldouble setTimeStep(bldouble CFL, bldouble dx, bldouble invVel) const;
+ bldouble computeSoundSpeed(bldouble gradU2) const;
+ void resetSolution(size_t inod, std::shared_ptr<BoundaryLayer> &bl, std::vector<bldouble> Sln0,
bool usePrevPoint = false);
};
} // namespace blast
diff --git a/blast/src/blast.h b/blast/src/blast.h
index 42b50f849775ecf7ade8d1e49a8690ca2bb3ed86..239c65fac8c8767324e5eed99021ce489cabce5a 100644
--- a/blast/src/blast.h
+++ b/blast/src/blast.h
@@ -30,6 +30,23 @@
#endif
#include "tbox.h"
+#include "blaster_def.h"
+#include <Eigen/Dense>
+
+#ifndef SWIG
+#ifndef USE_CODI
+using bldouble = double;
+#else // USE_CODI
+#include "codi.hpp"
+using bldouble = codi::RealReverse;
+using bltape = bldouble::Tape;
+#endif // USE_CODI
+
+using blMatrixXd = Eigen::Matrix<bldouble, Eigen::Dynamic, Eigen::Dynamic>;
+using blVectorXd = Eigen::Matrix<bldouble, Eigen::Dynamic, 1>;
+using blMatrix3d = Eigen::Matrix<bldouble, 3, 3>;
+using blVector3d = Eigen::Matrix<bldouble, 3, 1>;
+#endif // SWIG
/**
* @brief Namespace for blast module
diff --git a/blast/tests/dart/t_naca0012_2D.py b/blast/tests/dart/t_naca0012_2D.py
index b31c31d4d44a7699027636c45f7e59cde75ebd6a..2d169708ea564fa9e3a5ae1d2b2e4ffd53d8c1bc 100644
--- a/blast/tests/dart/t_naca0012_2D.py
+++ b/blast/tests/dart/t_naca0012_2D.py
@@ -121,7 +121,7 @@ def main():
# 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, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+ 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, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.getCdt(), vsol.getCdp(), vsol.getCdf(), isol.getCm()))
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
@@ -136,9 +136,9 @@ def main():
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('Cl', isol.getCl(), 0.230, 5e-2)) # XFOIL 0.2325
- tests.add(CTest('Cd wake', vsol.Cdt, 0.0057, 1e-3, forceabs=True)) # XFOIL 0.00531
- tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), 0.0058, 1e-3, forceabs=True)) # XFOIL 0.00531
- tests.add(CTest('Cdf', vsol.Cdf, 0.0047, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
+ tests.add(CTest('Cd wake', vsol.getCdt(), 0.0057, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cd integral', vsol.getCdf() + isol.getCd(), 0.0058, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cdf', vsol.getCdf(), 0.0047, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
tests.add(CTest('xtr_top', vsol.bodies[0].getAvgxtr(0), 0.196, 0.03, forceabs=True)) # XFOIL 0.1877
tests.add(CTest('xtr_bot', vsol.bodies[0].getAvgxtr(1), 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
tests.add(CTest('error dStar', np.linalg.norm(bl_solution['deltaStar'] - bl_solution['theta']*bl_solution['H']), 0., 0., forceabs=True))
diff --git a/blast/tests/dart/t_rae2822_3D.py b/blast/tests/dart/t_rae2822_3D.py
index b4018f567d9e44e42ed613a50b34ddf37826efa9..1bf4253032e9b757df58ef3d4f0d79f2e14b516a 100644
--- a/blast/tests/dart/t_rae2822_3D.py
+++ b/blast/tests/dart/t_rae2822_3D.py
@@ -143,7 +143,7 @@ def main():
# Display results.
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' Re M alpha Cl Cd Cd_wake 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} {8:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdf + isol.getCd(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+ 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} {8:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.getCdf() + isol.getCd(), vsol.getCdt(), vsol.getCdp(), vsol.getCdf(), isol.getCm()))
# Save results
isol.save(sfx='_viscous')
@@ -159,9 +159,9 @@ def main():
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('Cl', isol.getCl(), 0.1305, 5e-2))
- tests.add(CTest('Cd wake', vsol.Cdt, 0.0113, 1e-3, forceabs=True))
- tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), 0.0172, 1e-3, forceabs=True))
- tests.add(CTest('Cdf', vsol.Cdf, 0.0115, 1e-3, forceabs=True))
+ tests.add(CTest('Cd wake', vsol.getCdt(), 0.0113, 1e-3, forceabs=True))
+ tests.add(CTest('Cd integral', vsol.getCdf() + isol.getCd(), 0.0172, 1e-3, forceabs=True))
+ tests.add(CTest('Cdf', vsol.getCdf(), 0.0115, 1e-3, forceabs=True))
tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 4, 0, forceabs=True))
tests.run()
diff --git a/blaster_def.h.in b/blaster_def.h.in
new file mode 100644
index 0000000000000000000000000000000000000000..4d1fa8efc88a470e528b8a0146bcfc5b17ea07b6
--- /dev/null
+++ b/blaster_def.h.in
@@ -0,0 +1 @@
+#cmakedefine USE_CODI