diff --git a/dart/CMakeLists.txt b/dart/CMakeLists.txt
index b63c6eb6ba3e3bc90c0fdb2e56045d1aac07e87c..0ac4587bff42c5e93fce7ce74c10c5ead64ffe07 100644
--- a/dart/CMakeLists.txt
+++ b/dart/CMakeLists.txt
@@ -1,11 +1,11 @@
# Copyright 2020 University of Liège
-#
+#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
-#
+#
# http://www.apache.org/licenses/LICENSE-2.0
-#
+#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -23,6 +23,7 @@ MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
${CMAKE_CURRENT_LIST_DIR}/utils.py
DESTINATION dart)
-INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/scripts
+INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/interfaces
+ ${CMAKE_CURRENT_LIST_DIR}/scripts
${CMAKE_CURRENT_LIST_DIR}/viscous
DESTINATION dart)
diff --git a/dart/benchmark/onera.py b/dart/benchmark/onera.py
index 3c42b26029de2940079c66c0a8da79a0ec516caa..910b225a72cd1d677876f0e74301091afc28679c 100644
--- a/dart/benchmark/onera.py
+++ b/dart/benchmark/onera.py
@@ -40,21 +40,16 @@ def newton(pbl):
import dart
import tbox
# Pardiso and GMRES should give similar perfs
+ k = parseargs().k
try:
- newton = dart.Newton(pbl, tbox.Pardiso())
+ newton = dart.Newton(pbl, tbox.Pardiso(), nthrds=k, vrb=2)
except:
gmres = tbox.Gmres()
gmres.setFillFactor(2)
gmres.setDropTol(1e-6)
gmres.setRestart(50)
gmres.setTolerance(1e-5)
- newton = dart.Newton(pbl, gmres)
- newton.nthreads = parseargs().k
- newton.verbose = 1
- newton.relTol = 1e-6
- newton.absTol = 1e-8
- newton.lsTol = 1e-6
- newton.avThrsh = 1e-2
+ newton = dart.Newton(pbl, gmres, nthrds=k, vrb=2)
return newton
def main():
diff --git a/dart/default.py b/dart/default.py
index 6459935a68394e1e2913942bd007810a0ecf5d4f..e86e382a3b687495a0b0ce6e013cd21b90bae105 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -90,13 +90,7 @@ def picard(pbl):
'''Initialize Picard solver
'''
args = parseargs()
- solver = dart.Picard(pbl, tbox.Gmres())
- solver.nthreads = args.k
- solver.verbose = args.verb
- solver.relTol = 1e-6
- solver.absTol = 1e-8
- solver.maxIt = 25
- solver.relax = 0.7
+ solver = dart.Picard(pbl, tbox.Gmres(), nthrds=args.k, vrb=args.verb+1)
return solver
def newton(pbl):
@@ -104,15 +98,7 @@ def newton(pbl):
'''
from tbox.solvers import LinearSolver
args = parseargs()
- solver = dart.Newton(pbl, LinearSolver().pardiso())
- solver.nthreads = args.k
- solver.verbose = args.verb
- solver.relTol = 1e-6
- solver.absTol = 1e-8
- solver.maxIt = 25
- solver.lsTol = 1e-6
- solver.maxLsIt = 10
- solver.avThrsh = 1e-2
+ solver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb+1)
return solver
def meshMorpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = 'field', wk = 'wake', sym = 'symmetry'):
diff --git a/dart/interfaces/__init__.py b/dart/interfaces/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..40a96afc6ff09d58a702b76e3f7dd412fe975e26
--- /dev/null
+++ b/dart/interfaces/__init__.py
@@ -0,0 +1 @@
+# -*- coding: utf-8 -*-
diff --git a/dart/interfaces/mphys.py b/dart/interfaces/mphys.py
new file mode 100644
index 0000000000000000000000000000000000000000..64591d5aba34e11cf88bf2e789ca5f9c1239c372
--- /dev/null
+++ b/dart/interfaces/mphys.py
@@ -0,0 +1,647 @@
+# -*- coding: utf-8 -*-
+
+# Copyright 2021 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.
+
+## MPHYS interface
+# Adrien Crovato
+#
+# Structure:
+# - AeroBuilder
+# - AeroSurfaceMesh
+# - AeroGroup
+# - Morpher
+# - Solver
+# - Loads
+# - AeroCoefficients
+
+import numpy as np
+import openmdao.api as om
+from mphys.builder import Builder
+
+# Surface mesh
+class DartMesh(om.IndepVarComp):
+ """Initial surface mesh of moving boundary
+ For compatibility, the node coordinates of connected surfaces are always 3D
+
+ Attributes
+ ----------
+ dim : int
+ Problem dimension
+ bnd : dart.Boundary object
+ Moving boundary
+ """
+ def initialize(self):
+ self.options.declare('dim', desc='problem dimension')
+ self.options.declare('bnd', desc='moving boundary', recordable=False)
+
+ def setup(self):
+ self.bnd = self.options['bnd']
+ # Store all the node coordinates and add as output
+ x_aero0 = np.zeros(3 * len(self.bnd.nodes)) # surface node coordinates are 3D
+ for i in range(len(self.bnd.nodes)):
+ for j in range(3):
+ x_aero0[3 * i + j] = self.bnd.nodes[i].pos[j]
+ self.add_output('x_aero0', val=x_aero0, desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates'])
+
+ def mphys_get_triangulated_surface(self):
+ """Triangulate the surface
+ Only for 2D surfaces in 3D problems
+ Return a list containing a point and two vectors for each surface element
+ o p0
+ / \
+ v0 /___\ v1
+ """
+ if self.options['dim'] == 2:
+ raise RuntimeError('DartMesh - cannot triangulate 1D surface!\n')
+ # Create list of point coordinates and vector components
+ p0 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
+ v0 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
+ v1 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
+ for i in range(len(self.bnd.groups[0].tag.elems)):
+ e = self.bnd.groups[0].tag.elems[i]
+ p0i = e.nodes[0].pos
+ v0i = e.nodes[1].pos - e.nodes[0].pos
+ v1i = e.nodes[2].pos - e.nodes[0].pos
+ for j in range(3):
+ p0[i][j] = p0i[j]
+ v0[i][j] = v0i[j]
+ v1[i][j] = v1i[j]
+ return [p0, v0, v1]
+
+# Mesh morphers
+class DartMorpher(om.ImplicitComponent):
+ """Volume mesh morpher
+ For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
+
+ Attributes
+ ----------
+ dim : int
+ Problem dimension, used to set the length of the volume node coordinates vector
+ bnd : dart.Boundary object
+ Moving boundary
+ mrf : tbox.MshDeform object
+ Mesh morpher
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('dim', desc='problem dimension')
+ self.options.declare('bnd', desc='moving boundary', recordable=False)
+ self.options.declare('mrf', desc='mesh morpher', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ self.dim = self.options['dim']
+ self.bnd = self.options['bnd']
+ self.mrf = self.options['mrf']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
+ self.add_output('xv', val=np.zeros(self.dim * len(self.mrf.msh.nodes)), desc='aerodynamic volume node coordinates', tags=['mphys_coupling']) # volume node coordinates can be 2D or 3D
+ # Partials
+ # self.declare_partials(of=['xv'], wrt=['x_aero'])
+
+ def solve_nonlinear(self, inputs, outputs):
+ """Deform the volume mesh
+ """
+ # Update inputs
+ self.mrf.savePos()
+ for i in range(len(self.bnd.nodes)):
+ for j in range(3):
+ self.bnd.nodes[i].pos[j] = inputs['x_aero'][3 * i + j]
+ # Compute
+ self.mrf.deform()
+ # Update outputs
+ for i in range(len(self.mrf.msh.nodes)):
+ for j in range(self.dim):
+ outputs['xv'][self.dim * i + j] = self.mrf.msh.nodes[i].pos[j]
+
+ def apply_nonlinear(self, inputs, outputs, residuals):
+ """Compute the residuals of the mesh coordinates
+ """
+ raise NotImplementedError('DartMorpher - cannot compute the residuals!\n')
+
+ def linearize(self, inputs, outputs, partials):
+ """Compute the mesh jacobian matrix
+ """
+ self.adj.omLinearizeMesh()
+
+ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
+ """Compute the partials derivatives and perform the matrix-vector product
+ """
+ if mode == 'rev':
+ if 'xv' in d_residuals:
+ if 'xv' in d_outputs:
+ d_outputs['xv'] += self.adj.omComputeJacVecRmeshMesh(d_residuals['xv']) # dX = dRx_dX^T * dRx
+ if 'x_aero' in d_inputs:
+ for i in range(len(self.bnd.nodes)):
+ for j in range(self.dim):
+ d_inputs['x_aero'][3 * i + j] -= d_residuals['xv'][self.dim * self.bnd.nodes[i].row + j] # dXs = -dX
+ elif mode == 'fwd':
+ raise NotImplementedError('DartMorpher - forward mode not implemented!\n')
+
+ def solve_linear(self, d_outputs, d_residuals, mode):
+ """Solve for the partials gradients of the mesh residuals
+ """
+ if mode == 'rev':
+ d_residuals['xv'] = self.adj.omSolveMesh(d_outputs['xv']) # dRx = dRx_dX^-T * dX
+ elif mode == 'fwd':
+ raise NotImplementedError('DartMorpher - forward mode not implemented!\n')
+
+class DartDummyMorpher(om.ExplicitComponent):
+ """Dummy volume mesh morpher
+ Allow to define and link the surface mesh coordinates to the volume mesh coordinates, in case the mesh do not deform.
+ For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
+ """
+ def initialize(self):
+ self.options.declare('dim', desc='problem dimension')
+ self.options.declare('msh', desc='mesh data structure', recordable=False)
+
+ def setup(self):
+ dim = self.options['dim']
+ msh = self.options['msh']
+ # Set initial values
+ xv = np.zeros(dim * len(msh.nodes)) # volume node coordinates can be 2D or 3D
+ for i in range(len(msh.nodes)):
+ for j in range(dim):
+ xv[dim * i + j] = msh.nodes[i].pos[j]
+ # I/O
+ self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
+ self.add_output('xv', val=xv, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+
+ def compute(self, inputs, outputs):
+ """DO NOT deform the volume mesh
+ """
+ pass
+
+ def compute_partials(self, inputs, partials):
+ """DO NOT compute the partials derivatives
+ """
+ pass
+
+# Solver
+class DartSolver(om.ImplicitComponent):
+ """Aerodynamic solver
+
+ Attributes
+ ----------
+ sol : dart.Newton object
+ Direct Newton solver
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ self.sol = self.options['sol']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('aoa', val=0., units='rad', desc='angle of attack', tags=['mphys_input'])
+ self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+ self.add_output('phi', val=np.zeros(len(self.sol.pbl.msh.nodes)), desc='dart variables (potential)', tags=['mphys_coupling'])
+ # Partials
+ # self.declare_partials(of=['phi'], wrt=['aoa', 'xv', 'phi'])
+
+ def solve_nonlinear(self, inputs, outputs):
+ """Compute the flow variables
+ """
+ # Update inputs
+ # NB: the mesh is already up-to-date mesh already up-to-date because DartMorpher has already been run
+ self.sol.pbl.update(inputs['aoa'][0])
+ # Compute
+ self.sol.reset() # resets ICs (slows down convergence, but increases robustness for aero-structural)
+ self.sol.run()
+ # Update outputs
+ outputs['phi'] = self.sol.phi
+
+ def apply_nonlinear(self, inputs, outputs, residuals):
+ """Compute the residuals of the flow variables
+ """
+ raise NotImplementedError('DartSolver - cannot compute the residuals!\n')
+
+ def linearize(self, inputs, outputs, partials):
+ """Compute the flow jacobian matrix
+ """
+ self.adj.omLinearizeFlow()
+
+ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
+ """Compute the partials derivatives and perform the matrix-vector product
+ """
+ if mode == 'rev':
+ if 'phi' in d_residuals:
+ if 'phi' in d_outputs:
+ d_outputs['phi'] += self.adj.omComputeJacVecResidualFlow(d_residuals['phi']) # dPhi = dRphi_dPhi^T * dRphi
+ if 'aoa' in d_inputs:
+ d_inputs['aoa'] += self.adj.omComputeJacVecResidualAoa(d_residuals['phi']) # dAoa = dRphi_dAoa^T * dRphi
+ if 'xv' in d_inputs:
+ d_inputs['xv'] += self.adj.omComputeJacVecResidualMesh(d_residuals['phi']) # dX = dRphi_dX^T * dRphi
+ elif mode == 'fwd':
+ raise NotImplementedError('DartSolver - forward mode not implemented!\n')
+
+ def solve_linear(self, d_outputs, d_residuals, mode):
+ """Solve for the partials gradients of the flow residuals
+ """
+ if mode == 'rev':
+ d_residuals['phi'] = self.adj.omSolveFlow(d_outputs['phi']) # dRphi = dRphi_dPhi^-T * dPhi
+ elif mode == 'fwd':
+ raise NotImplementedError('DartSolver - forward mode not implemented!\n')
+
+# Aerodynamic loads
+class DartLoads(om.ExplicitComponent):
+ """Aerodynamic loads on moving boundary
+ For compatibility, the forces of connected surfaces are always 3D
+
+ Attributes
+ ----------
+ qinf : dart.Newton object
+ Freestream dynamic pressure
+ bnd : dart.Boundary object
+ Moving boundary
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('qinf', desc='freestream dynamic pressure')
+ self.options.declare('bnd', desc='moving boundary', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ self.qinf = self.options['qinf']
+ self.bnd = self.options['bnd']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+ self.add_input('phi', shape_by_conn=True, desc='flow variables (potential)', tags=['mphys_coupling'])
+ self.add_output('f_aero', val=np.zeros(3 * len(self.bnd.nodes)), desc='aerodynamic loads', tags=['mphys_coupling'])
+ # Partials
+ # self.declare_partials(of=['f_aero'], wrt=['xv', 'phi'])
+
+ def compute(self, inputs, outputs):
+ """Get the forces on moving boundary
+ """
+ # NB: inputs already up-to-date because DartSolver has already been run
+ for i in range(len(self.bnd.nodes)):
+ for j in range(3):
+ outputs['f_aero'][3 * i + j] = self.qinf * self.bnd.nLoads[i][j]
+
+ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
+ """Compute the partials derivatives and perform the matrix-vector product
+ """
+ if mode == 'rev':
+ if 'f_aero' in d_outputs:
+ if 'xv' in d_inputs:
+ d_inputs['xv'] += self.qinf * np.asarray(self.adj.omComputeJacVecLoadsMesh(d_outputs['f_aero'], self.bnd)) # dX = dL_dX^T * dL
+ if 'phi' in d_inputs:
+ d_inputs['phi'] += self.qinf * np.asarray(self.adj.omComputeJacVecLoadsFlow(d_outputs['f_aero'], self.bnd)) # dPhi = dL_dPhi^T * dL
+ elif mode == 'fwd':
+ raise NotImplementedError('DartLoads - forward mode not implemented!\n')
+
+# Aerodynamic coefficients
+class DartCoefficients(om.ExplicitComponent):
+ """Aerodynamic load coefficients for full aircraft
+ Also write solution to disk
+
+ Attributes
+ ----------
+ iscn : int
+ ID of scenario (default: 0)
+ sol : dart.Newton object
+ Direct Newton sovler
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+ self.options.declare('wrtr', desc='data writer', recordable=False)
+ self.options.declare('morph', default=False, desc='whether the mesh can deform or not')
+
+ def setup(self):
+ self.iscn = 0 # id of scenario
+ self.sol = self.options['sol']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('aoa', val=0., units='rad', desc='angle of attack', tags=['mphys_input'])
+ self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+ self.add_input('phi', shape_by_conn=True, desc='flow variables (potential)', tags=['mphys_coupling'])
+ self.add_output('cl', val=0., desc='lift coefficient', tags=['mphys_result'])
+ self.add_output('cd', val=0., desc='drag coefficient', tags=['mphys_result'])
+ # Partials
+ self.declare_partials(of=['cl', 'cd'], wrt=['aoa', 'xv', 'phi'])
+
+ def mphys_set_iscn(self, iscn):
+ """Set the id of the scenario
+ Allow to save the mesh and results in different files
+ """
+ self.iscn = iscn
+
+ def compute(self, inputs, outputs):
+ """Get the coefficients for the full aircraft
+ """
+ # NB: inputs already up-to-date because DartSolver has already been run
+ # Write to disk
+ if self.options['morph'] and str(type(self.options['wrtr'])) == '<class \'tboxw.GmshExport\'>':
+ if self.iscn == 0:
+ self.options['wrtr'].save(self.sol.pbl.msh.name)
+ else:
+ self.options['wrtr'].save(self.sol.pbl.msh.name + '_{0:d}'.format(self.iscn))
+ self.sol.save(self.options['wrtr'], self.iscn)
+ # Update outputs
+ outputs['cl'] = self.sol.Cl
+ outputs['cd'] = self.sol.Cd
+
+ def compute_partials(self, inputs, partials):
+ """Compute the partials derivatives
+ """
+ # Gradients wrt to angle of attack
+ dCfAoa = self.adj.omComputeGradientCoefficientsAoa()
+ partials['cl', 'aoa'] = dCfAoa[0]
+ partials['cd', 'aoa'] = dCfAoa[1]
+ # Gradients wrt to mesh
+ dCfMsh = self.adj.omComputeGradientCoefficientsMesh()
+ partials['cl', 'xv'] = dCfMsh[0]
+ partials['cd', 'xv'] = dCfMsh[1]
+ # Gradients wrt to flow variables
+ dCfPhi = self.adj.omComputeGradientCoefficientsFlow()
+ partials['cl', 'phi'] = dCfPhi[0]
+ partials['cd', 'phi'] = dCfPhi[1]
+
+# Aerodynamic group
+class DartGroup(om.Group):
+ """Aerodynamic group
+ Integrate the aerodynamic computations in an aero-structural coupling procedure
+
+ Components
+ ----------
+ - Mesh morpher
+ - Direct and adjoint solvers
+ - Loads computation
+ """
+ def initialize(self):
+ self.options.declare('qinf', desc='freestream dynamic pressure')
+ self.options.declare('bnd', desc='moving boundary', recordable=False)
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('mrf', desc='mesh morpher', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ # Components
+ if self.options['mrf'] is None:
+ self.add_subsystem('morpher', DartDummyMorpher(dim = self.options['sol'].pbl.nDim, msh = self.options['sol'].pbl.msh), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
+ else:
+ self.add_subsystem('morpher', DartMorpher(dim = self.options['sol'].pbl.nDim, bnd = self.options['bnd'], mrf = self.options['mrf'], adj = self.options['adj']), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
+ self.add_subsystem('solver', DartSolver(sol = self.options['sol'], adj = self.options['adj']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
+ if self.options['qinf'] != 0:
+ self.add_subsystem('loads', DartLoads(qinf = self.options['qinf'], bnd = self.options['bnd'], adj = self.options['adj']), promotes_inputs=['xv', 'phi'], promotes_outputs=['f_aero'])
+
+# Builder
+class DartBuilder(Builder):
+ """Dart builder for MPHYS
+
+ Attributes
+ ----------
+ __dim : int
+ Dimension of the problem
+ __qinf : float
+ Freestream dynamic pressure
+ __msh : tbox.MshData object
+ Mesh data structure
+ __wrtr : tbox.MshExport object
+ Mesh writer
+ __mrf : tbox.MshDeform object
+ Mesh morpher
+ __pbl : dart.Probem object
+ Problem definition
+ __sol : dart.Newton object
+ Direct Newton solver
+ __adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def __init__(self, params, scenario='aerodynamic', task='analysis'):
+ """Instantiate and initialize all the solver components.
+ Because we do not use MPI, we do not use the initialize method.
+
+ Parameters
+ ----------
+ params : dict
+ Dictonnary of parameters for the mesh, solver, ...
+ scenario : string, optional
+ Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
+ task : string, optional
+ Type of task (available: analysis, optimization) (default: analysis)
+ """
+ # Initialize
+ print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
+
+ # Imports
+ from dartflo import tbox
+ from dartflo.tbox import gmsh
+ from dartflo import dart
+
+ # Basic checks and config
+ # dimension
+ if params['Dim'] != 2 and params['Dim'] != 3:
+ raise RuntimeError('Problem dimension should be 2 or 3, but ' + params['Dim'] + ' was given!\n')
+ else:
+ self.__dim = params['Dim']
+ # aoa and aos
+ if 'AoA' in params:
+ alpha = params['AoA'] * np.pi / 180
+ else:
+ alpha = 0
+ if self.__dim == 2:
+ if 'AoS' in params and params['AoS'] != 0:
+ raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
+ else:
+ beta = 0
+ if 'Symmetry' in params:
+ raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
+ else:
+ if 'AoS' in params:
+ if params['AoS'] != 0 and 'Symmetry' in params:
+ raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
+ beta = params['AoS'] * np.pi / 180
+ else:
+ beta = 0
+ # save format
+ if params['Format'] == 'vtk':
+ try:
+ import tboxVtk
+ Writer = tboxVtk.VtkExport
+ print('VTK libraries found! Results will be saved in VTK format.\n')
+ except:
+ Writer = tbox.GmshExport
+ print('VTK libraries not found! Results will be saved in gmsh format.\n')
+ else:
+ Writer = tbox.GmshExport
+ print('Results will be saved in gmsh format.\n')
+ # number of threads
+ if 'Threads' in params:
+ nthrd = params['Threads']
+ else:
+ nthrd = 1
+ # verbosity
+ if 'Verb' in params:
+ verb = params['Verb']
+ else:
+ verb = 0
+ # scenario and task type
+ if scenario != 'aerodynamic' and scenario != 'aerostructural':
+ raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n')
+ if task != 'analysis' and task != 'optimization':
+ raise RuntimeError('Task should be analysis or optimization, but "' + scenario + '" was given!\n')
+ # dynamic pressure
+ if scenario == 'aerostructural':
+ self.__qinf = params['Qinf']
+ else:
+ self.__qinf = 0
+
+ # Mesh creation
+ self.__msh = gmsh.MeshLoader(params['File']).execute(**params['Pars'])
+ # add the wake
+ if self.__dim == 2:
+ mshCrck = tbox.MshCrack(self.__msh, self.__dim)
+ mshCrck.setCrack(params['Wake'])
+ mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wing']])
+ mshCrck.run()
+ else:
+ for i in range(len(params['Wakes'])):
+ mshCrck = tbox.MshCrack(self.__msh, self.__dim)
+ mshCrck.setCrack(params['Wakes'][i])
+ mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wings'][i]])
+ if 'Fuselage' in params:
+ mshCrck.addBoundaries([params['Fuselage']])
+ if 'Symmetry' in params:
+ mshCrck.addBoundaries([params['Symmetry']])
+ if 'WakeTips' in params:
+ mshCrck.setExcluded(params['WakeTips'][i]) # 2.5D computations
+ mshCrck.run()
+ # save the updated mesh in native (gmsh) format and then set the writer to the requested format
+ tbox.GmshExport(self.__msh).save(self.__msh.name)
+ del mshCrck
+ self.__wrtr = Writer(self.__msh)
+ print('\n')
+
+ # Mesh morpher creation
+ if scenario == 'aerostructural' or task == 'optimization':
+ self.__mrf = tbox.MshDeform(self.__msh, self.__dim)
+ self.__mrf.nthreads = nthrd
+ self.__mrf.setField(params['Fluid'])
+ self.__mrf.addFixed(params['Farfield'])
+ if self.__dim == 2:
+ self.__mrf.addMoving([params['Wing']])
+ self.__mrf.addInternal([params['Wake'], params['Wake']+'_'])
+ else:
+ if 'Fuselage' in params:
+ self.__mrf.addFixed(params['Fuselage'])
+ self.__mrf.setSymmetry(params['Symmetry'], 1)
+ for i in range(len(params['Wings'])):
+ if i == 0:
+ self.__mrf.addMoving([params['Wings'][i]]) # TODO boundary of interest (FSI/OPT) is always first given boundary
+ else:
+ self.__mrf.addFixed([params['Wings'][i]])
+ self.__mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
+ self.__mrf.initialize()
+ else:
+ self.__mrf = None
+
+ # Problem creation
+ self.__pbl = dart.Problem(self.__msh, self.__dim, alpha, beta, params['M_inf'], params['S_ref'], params['c_ref'], params['x_ref'], params['y_ref'], params['z_ref'])
+ # add the field
+ self.__pbl.set(dart.Medium(self.__msh, params['Fluid'], params['M_inf'], self.__dim, alpha, beta))
+ # add the initial condition
+ self.__pbl.set(dart.Initial(self.__msh, params['Fluid'], self.__dim, alpha, beta))
+ # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
+ for i in range (len(params['Farfield'])):
+ self.__pbl.add(dart.Freestream(self.__msh, params['Farfield'][i], self.__dim, alpha, beta))
+ # add solid boundaries
+ if self.__dim == 2:
+ bnd = dart.Boundary(self.__msh, [params['Wing'], params['Fluid']])
+ self.__bnd = bnd
+ self.__pbl.add(bnd)
+ else:
+ self.__bnd = None
+ for bd in params['Wings']:
+ bnd = dart.Boundary(self.__msh, [bd, params['Fluid']])
+ if self.__bnd is None:
+ self.__bnd = bnd # TODO boundary of interest (FSI/OPT) is always first given boundary
+ self.__pbl.add(bnd)
+ if 'Fuselage' in params:
+ bnd = dart.Boundary(self.__msh, [params['Fuselage'], params['Fluid']])
+ self.__pbl.add(bnd)
+ # add Wake/Kutta boundary conditions
+ if self.__dim == 2:
+ self.__pbl.add(dart.Wake(self.__msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
+ self.__pbl.add(dart.Kutta(self.__msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
+ else:
+ for i in range(len(params['Wakes'])):
+ self.__pbl.add(dart.Wake(self.__msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['TeTips'][i]]))
+
+ # Direct (forward) solver creation
+ # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones
+ # initialize the linear (inner) solver
+ if params['LSolver'] == 'PARDISO':
+ linsol = tbox.Pardiso()
+ elif params['LSolver'] == 'MUMPS':
+ linsol = tbox.Mumps()
+ elif params['LSolver'] == 'GMRES':
+ linsol = tbox.Gmres()
+ linsol.setFillFactor(params['G_fill']) if 'G_fill' in params else linsol.setFillFactor(2)
+ linsol.setRestart(params['G_restart']) if 'G_restart' in params else linsol.setRestart(50)
+ linsol.setTolerance(params['G_tol']) if 'G_tol' in params else linsol.setTolerance(1e-5)
+ else:
+ raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + params['LSolver'] + ' was given!\n')
+ # initialize the nonlinear (outer) solver
+ self.__sol = dart.Newton(self.__pbl, linsol, rTol=params['Rel_tol'], aTol=params['Abs_tol'], mIt=params['Max_it'], nthrds=nthrd, vrb=verb)
+
+ # Adjoint (reverse) solver creation
+ if task == 'optimization':
+ if params['LSolver'] == 'PARDISO':
+ alinsol = tbox.Pardiso()
+ elif params['LSolver'] == 'MUMPS':
+ alinsol = tbox.Mumps()
+ else:
+ alinsol = tbox.Gmres()
+ alinsol.setFillFactor(params['G_fill']) if 'G_fill' in params else alinsol.setFillFactor(2)
+ alinsol.setRestart(params['G_restart']) if 'G_restart' in params else alinsol.setRestart(50)
+ alinsol.setTolerance(1e-12)
+ self.__adj = dart.Adjoint(self.__sol, self.__mrf, alinsol)
+ else:
+ self.__adj = None
+
+ # Finalize
+ print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
+
+ def get_mesh_coordinate_subsystem(self):
+ """Return openMDAO component to get the initial surface mesh coordinates
+ """
+ return DartMesh(dim = self.__dim, bnd = self.__bnd)
+
+ def get_coupling_group_subsystem(self):
+ """Return openMDAO group containing the morpher, solver and loads
+ """
+ return DartGroup(qinf = self.__qinf, bnd = self.__bnd, sol = self.__sol, mrf = self.__mrf, adj = self.__adj)
+
+ def get_post_coupling_subsystem(self):
+ """Return openMDAO component that computes the aero coefficients and writes data to disk
+ """
+ return DartCoefficients(sol = self.__sol, adj = self.__adj, wrtr = self.__wrtr, morph = False if self.__mrf is None else True)
+
+ def get_number_of_nodes(self):
+ """Return the number of surface nodes
+ """
+ return len(self.__bnd.nodes)
diff --git a/dart/scripts/config.py b/dart/scripts/config.py
index 00b934863d7c815482e0e3ec56d4a4d84df01915..4ba055dbb1ebeff0350d2e9b53adcd5f9bf0ef7e 100644
--- a/dart/scripts/config.py
+++ b/dart/scripts/config.py
@@ -134,19 +134,14 @@ class Config:
else:
raise RuntimeError('Available linear solvers: Pardiso, GMRES, MUMPS or SparseLU, but ' + p['LSolver'] + ' was given!\n')
# initialize the nonlinear (outer) solver
+ args = parseargs()
if p['NSolver'] == 'Picard':
- self.solver = dart.Picard(self.pbl, linsol)
+ self.solver = dart.Picard(self.pbl, linsol, rTol=p['Rel_tol'], aTol=p['Abs_tol'], mIt=p['Max_it'], nthrds=args.k, vrb=args.verb+1)
self.solver.relax = p['Relaxation']
elif p['NSolver'] == 'Newton':
- self.solver = dart.Newton(self.pbl, linsol)
+ self.solver = dart.Newton(self.pbl, linsol, rTol=p['Rel_tol'], aTol=p['Abs_tol'], mIt=p['Max_it'], nthrds=args.k, vrb=args.verb+1)
self.solver.lsTol = p['LS_tol']
self.solver.maxLsIt = p['Max_it_LS']
self.solver.avThrsh = p['AV_thrsh']
else:
raise RuntimeError('Available nonlinear solvers: Picard or Newton, but ' + p['NSolver'] + ' was given!\n')
- args = parseargs()
- self.solver.nthreads = args.k
- self.solver.verbose = args.verb
- self.solver.relTol = p['Rel_tol']
- self.solver.absTol = p['Abs_tol']
- self.solver.maxIt = p['Max_it']
diff --git a/dart/src/wNewton.cpp b/dart/src/wNewton.cpp
index caa9e016a152a93b1b612fcbd1530c1c24ee79a6..ce523f85741754783d2a699dfb2e68b70a3b1e13 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -57,8 +57,18 @@ using namespace dart;
* @brief Initialize the solver and set default parameters
* @authors Adrien Crovato
*/
-Newton::Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol) : Solver(_pbl, _linsol)
+Newton::Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol,
+ double rTol, double aTol, int mIt,
+ int nthrds, int vrb) : Solver(_pbl, _linsol, rTol, aTol, mIt, nthrds, vrb)
{
+ // Additional display
+ std::cout << "--- Newton Solver ---\n"
+ << "Inner solver: " << *linsol
+ << "Number of threads: " << nthreads << "\n"
+ << "Relative tolerance: " << log10(relTol) << "\n"
+ << "Absolute tolerance: " << log10(absTol) << "\n"
+ << "Maximum iterations: " << maxIt << "\n"
+ << std::endl;
// Additional default parameters
lsTol = 1e-6;
maxLsIt = 10;
@@ -78,23 +88,12 @@ bool Newton::run()
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
- // Display solver and problem parameters
- std::cout << "--- Newton solver ---\n"
- << "Inner solver: " << *linsol
- << std::setprecision(5)
- << "Number of threads: " << nthreads << "\n"
- << "Relative tolerance: " << log10(relTol) << "\n"
- << "Absolute tolerance: " << log10(absTol) << "\n"
- << "Maximum iterations: " << maxIt << "\n";
- std::cout << "--- Problem definition ---\n"
- << std::setprecision(5)
- << "Angle of attack: " << pbl->alpha * 180 / 3.14159 << " deg\n"
- << "Angle of sideslip: " << pbl->beta * 180 / 3.14159 << " deg\n"
- << "Mach number: " << pbl->M_inf << "\n"
- << "Reference area: " << pbl->S_ref << "\n"
- << "Reference length: " << pbl->c_ref << "\n"
- << "Reference origin: (" << pbl->x_ref(0) << ", " << pbl->x_ref(1) << ", " << pbl->x_ref(2) << ")\n"
- << std::endl;
+ // Display current freestream conditions
+ std:: cout << std::setprecision(2)
+ << "- Mach " << pbl->M_inf << ", "
+ << pbl->alpha * 180 / 3.14159 << "° AoA, "
+ << pbl->beta * 180 / 3.14159 << "° AoS"
+ << std::endl;
// Initialize solver loop
nIt = 0; // iteration counter
@@ -199,15 +198,18 @@ bool Newton::run()
// Compute aerodynamic coefficient on boundaries
Solver::computeLoad();
- // Display residual
- std::cout << std::fixed << std::setprecision(5);
- std::cout << std::setw(8) << nIt
- << std::setw(8) << linsol->getIterations()
- << std::setw(8) << ls.fevalIt
- << std::setw(12) << Cl
- << std::setw(12) << Cd
- << std::setw(15) << log10(relRes)
- << std::setw(15) << log10(absRes) << "\n";
+ // Display residual (at each iteration)
+ if (verbose > 0)
+ {
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << nIt
+ << std::setw(8) << linsol->getIterations()
+ << std::setw(8) << ls.fevalIt
+ << std::setw(12) << Cl
+ << std::setw(12) << Cd
+ << std::setw(15) << log10(relRes)
+ << std::setw(15) << log10(absRes) << "\n";
+ }
// Check convergence
if ((relRes < relTol || absRes < absTol) && nIt != 0)
@@ -218,6 +220,19 @@ bool Newton::run()
continue;
} while (nIt < maxIt);
+ // Display residual (only last iteration)
+ if (verbose == 0)
+ {
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << nIt
+ << std::setw(8) << linsol->getIterations()
+ << std::setw(8) << ls.fevalIt
+ << std::setw(12) << Cl
+ << std::setw(12) << Cd
+ << std::setw(15) << log10(relRes)
+ << std::setw(15) << log10(absRes) << "\n";
+ }
+
// Compute field variables
Solver::computeFlow();
@@ -225,7 +240,7 @@ bool Newton::run()
if (relRes < relTol || absRes < absTol)
{
std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 0)
+ if (verbose > 1)
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
@@ -234,7 +249,7 @@ bool Newton::run()
else if (std::isnan(relRes))
{
std::cout << ANSI_COLOR_RED << "Newton solver diverged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 0)
+ if (verbose > 1)
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
@@ -243,7 +258,7 @@ bool Newton::run()
else
{
std::cout << ANSI_COLOR_YELLOW << "Newton solver not fully converged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 0)
+ if (verbose > 1)
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
@@ -378,7 +393,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
J.prune(0.);
J.makeCompressed();
- if (verbose > 1)
+ if (verbose > 2)
std::cout << "J (" << J.rows() << "," << J.cols() << ") nnz=" << J.nonZeros() << "\n";
}
@@ -487,7 +502,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
for (auto nod : dBC->nodes)
R(nod->row) = 0.;
- if (verbose > 1)
+ if (verbose > 2)
std::cout << "R (" << R.size() << ")\n";
}
diff --git a/dart/src/wNewton.h b/dart/src/wNewton.h
index 1a04e0cc1a07dcba7a8f616defb7d4000c4115b5..2499fac4102da25609ba2a3f2d1c67da20c83040 100644
--- a/dart/src/wNewton.h
+++ b/dart/src/wNewton.h
@@ -47,7 +47,7 @@ private:
double muCv; ///< variable artificial viscosity scaling factor
public:
- Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol);
+ Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol = 1e-6, double aTol = 1e-8, int mIt = 25, int nthrds = 1, int vrb = 1);
virtual bool run() override;
diff --git a/dart/src/wPicard.cpp b/dart/src/wPicard.cpp
index 9d38815fc5edd291d2f6509701eed1c8c3b0e93b..d503cec5529dd7c11e7099f940475c849ed9870a 100644
--- a/dart/src/wPicard.cpp
+++ b/dart/src/wPicard.cpp
@@ -55,10 +55,20 @@ using namespace dart;
* @brief Initialize Picard solver
* @authors Adrien Crovato
*/
-Picard::Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol) : Solver(_pbl, _linsol)
+Picard::Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol,
+ double rTol, double aTol, int mIt,
+ int nthrds, int vrb) : Solver(_pbl, _linsol, rTol, aTol, mIt, nthrds, vrb)
{
+ // Additional display
+ std::cout << "--- Picard Solver ---\n"
+ << "Inner solver: " << *linsol
+ << "Number of threads: " << nthreads << "\n"
+ << "Relative tolerance: " << log10(relTol) << "\n"
+ << "Absolute tolerance: " << log10(absTol) << "\n"
+ << "Maximum iterations: " << maxIt << "\n"
+ << std::endl;
// Additional default parameters
- relax = 1.;
+ relax = 0.7;
}
/**
@@ -72,23 +82,12 @@ bool Picard::run()
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
- // Display solver and problem parameters
- std::cout << "--- Picard solver ---\n"
- << "Inner solver: " << *linsol
- << std::setprecision(5)
- << "Number of threads: " << nthreads << "\n"
- << "Relative tolerance: " << log10(relTol) << "\n"
- << "Absolute tolerance: " << log10(absTol) << "\n"
- << "Maximum iterations: " << maxIt << "\n";
- std::cout << "--- Problem definition ---\n"
- << std::setprecision(5)
- << "Angle of attack: " << pbl->alpha * 180 / 3.14159 << " deg\n"
- << "Angle of sideslip: " << pbl->beta * 180 / 3.14159 << " deg\n"
- << "Mach number: " << pbl->M_inf << "\n"
- << "Reference area: " << pbl->S_ref << "\n"
- << "Reference length: " << pbl->c_ref << "\n"
- << "Reference origin: (" << pbl->x_ref(0) << ", " << pbl->x_ref(1) << ", " << pbl->x_ref(2) << ")\n"
- << std::endl;
+ // Display current freestream conditions
+ std:: cout << std::setprecision(2)
+ << "- Mach " << pbl->M_inf << ", "
+ << pbl->alpha * 180 / 3.14159 << "° AoA, "
+ << pbl->beta * 180 / 3.14159 << "° AoS"
+ << std::endl;
// Initialize solver loop
nIt = 0;
@@ -124,13 +123,16 @@ bool Picard::run()
// Compute aerodynamic coefficient on boundaries
Solver::computeLoad();
- // Display residual
- std::cout << std::fixed << std::setprecision(5);
- std::cout << std::setw(8) << nIt
- << std::setw(12) << Cl
- << std::setw(12) << Cd
- << std::setw(15) << log10(relRes)
- << std::setw(15) << log10(absRes) << "\n";
+ // Display residual (at each iteration)
+ if (verbose > 0 || nIt == 0)
+ {
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << nIt
+ << std::setw(12) << Cl
+ << std::setw(12) << Cd
+ << std::setw(15) << log10(relRes)
+ << std::setw(15) << log10(absRes) << "\n";
+ }
// Check convergence
if ((relRes < relTol || absRes < absTol) && nIt != 0)
@@ -151,6 +153,17 @@ bool Picard::run()
phi_ = (1 - relax) * phiOld + relax * phi_;
} while (nIt < maxIt);
+ // Display residual (only last iteration)
+ if (verbose == 0)
+ {
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << nIt
+ << std::setw(12) << Cl
+ << std::setw(12) << Cd
+ << std::setw(15) << log10(relRes)
+ << std::setw(15) << log10(absRes) << "\n";
+ }
+
// Compute field variables
Solver::computeFlow();
@@ -158,7 +171,7 @@ bool Picard::run()
if (relRes < relTol || absRes < absTol)
{
std::cout << ANSI_COLOR_GREEN << "Picard solver converged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 0)
+ if (verbose > 1)
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
@@ -167,7 +180,7 @@ bool Picard::run()
else if (std::isnan(relRes))
{
std::cout << ANSI_COLOR_RED << "Picard sovler diverged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 0)
+ if (verbose > 1)
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
@@ -176,7 +189,7 @@ bool Picard::run()
else
{
std::cout << ANSI_COLOR_YELLOW << "Picard solver not fully converged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 0)
+ if (verbose > 1)
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
@@ -332,7 +345,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
A.prune(0.);
A.makeCompressed();
- if (verbose > 1)
+ if (verbose > 2)
{
std::cout << "A (" << A.rows() << "," << A.cols() << ") nnz=" << A.nonZeros() << "\n";
std::cout << "b (" << b.size() << ")\n";
diff --git a/dart/src/wPicard.h b/dart/src/wPicard.h
index 7ff30338692e3c10bee9afb0e122578d75f4b8f9..3ebd9b88699479131a3ec4f2659adce9645bd116 100644
--- a/dart/src/wPicard.h
+++ b/dart/src/wPicard.h
@@ -38,7 +38,7 @@ public:
double relax; ///< relaxation
public:
- Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol);
+ Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol = 1e-6, double aTol = 1e-8, int mIt = 25, int nthrds = 1, int vrb = 1);
virtual bool run() override;
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index 1fc5bb3ed1bd1809a6b3b4a071233046875bafa8..78c87a44d3e414d21643711fe97550fbec864d5b 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -17,6 +17,7 @@
#include "wSolver.h"
#include "wProblem.h"
#include "wMedium.h"
+#include "wFreestream.h"
#include "wAssign.h"
#include "wBoundary.h"
#include "wWake.h"
@@ -26,12 +27,15 @@
#include "wNode.h"
#include "wElement.h"
#include "wTag.h"
+#include "wLinearSolver.h"
#include "wResults.h"
#include "wMshExport.h"
#include <tbb/parallel_for_each.h>
#include <tbb/spin_mutex.h>
+#include <iomanip>
+
using namespace tbox;
using namespace dart;
@@ -42,32 +46,30 @@ using namespace dart;
* @brief Initialize the solver and perform sanity checks
* @authors Adrien Crovato
*/
-Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _linsol) : pbl(_pbl), linsol(_linsol), nIt(0), Cl(0), Cd(0), Cs(0), Cm(0)
+Solver::Solver(std::shared_ptr<Problem> _pbl,
+ std::shared_ptr<LinearSolver> _linsol,
+ double rTol, double aTol, int mIt,
+ int nthrds, int vrb) : pbl(_pbl), linsol(_linsol),
+ relTol(rTol), absTol(aTol), maxIt(mIt),
+ nthreads(nthrds), verbose(vrb),
+ nIt(0), Cl(0), Cd(0), Cs(0), Cm(0)
{
// Say hi
- std::cout << std::endl;
- std::cout << "*******************************************************************************" << std::endl;
- std::cout << "** \\_/ **" << std::endl;
- std::cout << "** \\_______O(_)O_______/ **" << std::endl;
- std::cout << "** _______________________________ **" << std::endl;
- std::cout << "** ___ __ \\__ |__ __ \\__ __/ **" << std::endl;
- std::cout << "** __ / / /_ /| |_ /_/ /_ / **" << std::endl;
- std::cout << "** _ /_/ /_ ___ | _, _/_ / **" << std::endl;
- std::cout << "** /_____/ /_/ |_/_/ |_| /_/ **" << std::endl;
- std::cout << "*******************************************************************************" << std::endl;
- std::cout << "** Hi! My name is DART v1.0.0-21.09 **" << std::endl;
- std::cout << "** Adrien Crovato **" << std::endl;
- std::cout << "** ULiege 2018-2021 **" << std::endl;
- std::cout << "*******************************************************************************" << std::endl
+ std::cout << "*******************************************************************************\n"
+ << "** \\_/ **\n"
+ << "** \\_______O(_)O_______/ **\n"
+ << "** _______________________________ **\n"
+ << "** ___ __ \\__ |__ __ \\__ __/ **\n"
+ << "** __ / / /_ /| |_ /_/ /_ / **\n"
+ << "** _ /_/ /_ ___ | __ _/_ / **\n"
+ << "** /_____/ /_/ |_/_/ |_| /_/ **\n"
+ << "*******************************************************************************\n"
+ << "** Hi! My name is DART v1.1.0-21.09 **\n"
+ << "** Adrien Crovato **\n"
+ << "** ULiege 2018-2021 **\n"
+ << "*******************************************************************************\n"
<< std::endl;
- // Default parameters
- nthreads = 1;
- verbose = 0;
- maxIt = 50;
- relTol = 1e-6;
- absTol = 1e-8;
-
// Check the problem, pin a degree of freedom if necessary, and update element memory
pbl->check();
pbl->pin();
@@ -88,9 +90,7 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _lin
Cp.resize(pbl->msh->nodes.size(), 0.);
// Apply Initial and Dirichlet boundary conditions
- pbl->iIC->apply(phi);
- for (auto dBC : pbl->dBCs)
- dBC->apply(phi);
+ this->reset();
// Kutta condition step 1: equality of normal flux through the wake
// -> directly assemble upper wake contributions on lower wake rows, and leave the upper wake rows to zero
@@ -98,17 +98,54 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _lin
for (auto wake : pbl->wBCs)
for (auto p : wake->nodMap)
rows[p.first->row] = p.second->row;
+
+ // Display configuration
+ std::cout << "--- Mesh data ---\n"
+ << "Number of nodes: " << pbl->msh->nodes.size() << "\n"
+ << "Number of elements: " << pbl->msh->elems.size() << "\n";
+ std::cout << "--- Physical groups ---\n"
+ << "Fluid on: " << *pbl->medium->tag << "\n";
+ for (auto bnd : pbl->fBCs)
+ std::cout << "Freestream on " << *bnd->tag << "\n";
+ for (auto bnd : pbl->dBCs)
+ std::cout << "Dirichlet on " << *bnd->tag << "\n";
+ for (auto bnd : pbl->wBCs)
+ std::cout << "Wake on " << *bnd->groups[0]->tag << "\n";
+ for (auto bnd : pbl->bnds)
+ std::cout << "Boundary on " << *bnd->groups[0]->tag << "\n";
+ std::cout << "--- Freestream conditions ---\n"
+ << std::setprecision(5)
+ << "Angle of attack: " << pbl->alpha * 180 / 3.14159 << " deg\n"
+ << "Angle of sideslip: " << pbl->beta * 180 / 3.14159 << " deg\n"
+ << "Mach number: " << pbl->M_inf << "\n";
+ std::cout << "--- Reference data ---\n"
+ << "Reference area: " << pbl->S_ref << "\n"
+ << "Reference length: " << pbl->c_ref << "\n"
+ << "Reference origin: (" << pbl->x_ref(0) << ", " << pbl->x_ref(1) << ", " << pbl->x_ref(2) << ")"
+ << std::endl;
}
Solver::~Solver()
{
// Say bye
- std::cout << std::endl;
- std::cout << "*******************************************************************************" << std::endl;
- std::cout << "** \\_/ **" << std::endl;
- std::cout << "** \\_______O(_)O_______/ **" << std::endl;
- std::cout << "** ! ! **" << std::endl;
- std::cout << "*******************************************************************************" << std::endl;
- std::cout << std::endl;
+ std::cout << "\n"
+ << "*******************************************************************************\n"
+ << "** \\_/ **\n"
+ << "** \\_______O(_)O_______/ **\n"
+ << "** ! ! **\n"
+ << "*******************************************************************************\n"
+ << std::endl;
+}
+
+/**
+ * @brief Set or reset initial conditions
+ */
+void Solver::reset()
+{
+ // ICs
+ pbl->iIC->apply(phi);
+ // BCs
+ for (auto dBC : pbl->dBCs)
+ dBC->apply(phi);
}
/**
@@ -150,6 +187,7 @@ void Solver::save(MshExport *mshWriter, int n)
for (auto sur : pbl->bnds)
sur->save(sur->groups[0]->tag->name, results);
}
+ std::cout << std::endl;
}
/**
diff --git a/dart/src/wSolver.h b/dart/src/wSolver.h
index 0ccac4a701bf13f3b6352cb6287b3dc2d2af349f..f04d1bd5cfd74399e6058ec7614cb727897386d0 100644
--- a/dart/src/wSolver.h
+++ b/dart/src/wSolver.h
@@ -39,10 +39,10 @@ public:
std::shared_ptr<Problem> pbl; ///< problem definition
std::shared_ptr<tbox::LinearSolver> linsol; ///< linear solver
- int nthreads; ///< number of threads for the execution
double relTol; ///< relative tolerance on the residual
double absTol; ///< absolute tolerance on the residual
int maxIt; ///< max number of iterations
+ int nthreads; ///< number of threads for the execution
int verbose; ///< display more info
int nIt; ///< number of iterations
@@ -63,9 +63,10 @@ protected:
std::vector<int> rows; ///< unknown nodal index
public:
- Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol);
+ Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol, double aTol, int mIt, int nthrds, int vrb);
virtual ~Solver();
+ void reset();
virtual bool run();
void save(tbox::MshExport *mshWriter, int n = 0);