From 75ea7440a85265d03b414ad7da450e33deaf7768 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Mon, 20 Sep 2021 18:45:21 +0200
Subject: [PATCH] Common API (internal, cupydo, mphys)
---
dart/api/core.py | 233 ++++++++++++++++++++++++++++++++++++
dart/api/internal/config.py | 172 --------------------------
dart/api/internal/polar.py | 29 +++--
dart/api/internal/trim.py | 38 +++---
dart/api/mphys.py | 187 +----------------------------
dart/cases/coyote.py | 11 +-
dart/cases/n0012.py | 9 +-
dart/cases/n0012_3.py | 11 +-
dart/cases/n64A410.py | 11 +-
dart/cases/rae2822.py | 11 +-
dart/cases/wbht.py | 11 +-
dart/cases/wht.py | 11 +-
12 files changed, 318 insertions(+), 416 deletions(-)
create mode 100644 dart/api/core.py
delete mode 100644 dart/api/internal/config.py
diff --git a/dart/api/core.py b/dart/api/core.py
new file mode 100644
index 0000000..64cc2cf
--- /dev/null
+++ b/dart/api/core.py
@@ -0,0 +1,233 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Initialize DART
+# Adrien Crovato
+
+def initDart(params, scenario='aerodynamic', task='analysis'):
+ """Initlialize DART for various API
+
+ Parameters
+ ----------
+ params : dict
+ Dictionary of parameters to configure DART
+ scenario : string, optional
+ Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
+ task : string, optional
+ Type of task (available: analysis, optimization) (default: analysis)
+
+ Returns
+ -------
+ _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
+ _bnd : Dart.Body object
+ Body of interest
+ _sol : dart.Newton object
+ Direct (Newton) solver
+ _adj : dart.Adjoint object
+ Adjoint solver
+ """
+ # Imports
+ import math
+ import tbox
+ import tbox.gmsh as gmsh
+ 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:
+ _dim = params['Dim']
+ # aoa and aos
+ if 'AoA' in params:
+ alpha = params['AoA'] * math.pi / 180
+ else:
+ alpha = 0
+ if _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'] * math.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 "' + task + '" was given!\n')
+ # dynamic pressure
+ if scenario == 'aerostructural':
+ _qinf = params['Q_inf']
+ else:
+ _qinf = 0
+
+ # Mesh creation
+ _msh = gmsh.MeshLoader(params['File']).execute(**params['Pars'])
+ # add the wake
+ if _dim == 2:
+ mshCrck = tbox.MshCrack(_msh, _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(_msh, _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(_msh).save(_msh.name)
+ del mshCrck
+ _wrtr = Writer(_msh)
+ print('\n')
+
+ # Mesh morpher creation
+ if scenario == 'aerostructural' or task == 'optimization':
+ _mrf = tbox.MshDeform(_msh, _dim)
+ _mrf.nthreads = nthrd
+ _mrf.setField(params['Fluid'])
+ _mrf.addFixed(params['Farfield'])
+ if _dim == 2:
+ _mrf.addMoving([params['Wing']])
+ _mrf.addInternal([params['Wake'], params['Wake']+'_'])
+ else:
+ if 'Fuselage' in params:
+ _mrf.addFixed(params['Fuselage'])
+ _mrf.setSymmetry(params['Symmetry'], 1)
+ for i in range(len(params['Wings'])):
+ if i == 0:
+ _mrf.addMoving([params['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
+ else:
+ _mrf.addFixed([params['Wings'][i]])
+ _mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
+ _mrf.initialize()
+ else:
+ _mrf = None
+
+ # Problem creation
+ _pbl = dart.Problem(_msh, _dim, alpha, beta, params['M_inf'], params['S_ref'], params['c_ref'], params['x_ref'], params['y_ref'], params['z_ref'])
+ # add the field
+ _pbl.set(dart.Fluid(_msh, params['Fluid'], params['M_inf'], _dim, alpha, beta))
+ # add the initial condition
+ _pbl.set(dart.Initial(_msh, params['Fluid'], _dim, alpha, beta))
+ # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
+ for i in range (len(params['Farfield'])):
+ _pbl.add(dart.Freestream(_msh, params['Farfield'][i], _dim, alpha, beta))
+ # add solid boundaries
+ if _dim == 2:
+ bnd = dart.Body(_msh, [params['Wing'], params['Fluid']])
+ _bnd = bnd
+ _pbl.add(bnd)
+ else:
+ _bnd = None
+ for bd in params['Wings']:
+ bnd = dart.Body(_msh, [bd, params['Fluid']])
+ if _bnd is None:
+ _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
+ _pbl.add(bnd)
+ if 'Fuselage' in params:
+ bnd = dart.Body(_msh, [params['Fuselage'], params['Fluid']])
+ _pbl.add(bnd)
+ # add Wake/Kutta boundary conditions
+ if _dim == 2:
+ _pbl.add(dart.Wake(_msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
+ _pbl.add(dart.Kutta(_msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
+ else:
+ for i in range(len(params['Wakes'])):
+ _pbl.add(dart.Wake(_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
+ _sol = dart.Newton(_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)
+ _adj = dart.Adjoint(_sol, _mrf, alinsol)
+ else:
+ _adj = None
+
+ # Return
+ return _dim, _qinf, _msh, _wrtr, _mrf, _pbl, _bnd, _sol, _adj
diff --git a/dart/api/internal/config.py b/dart/api/internal/config.py
deleted file mode 100644
index 7c40318..0000000
--- a/dart/api/internal/config.py
+++ /dev/null
@@ -1,172 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## Base class for configuration
-# Adrien Crovato
-
-import math
-import dart
-import tbox
-import tbox.gmsh as gmsh
-import fwk
-from tbox.solvers import LinearSolver
-from fwk.wutils import parseargs
-from fwk.coloring import ccolors
-
-class Config:
- """Basic configurator
- Configure DART to perform standard computations on standard configurations
-
- Attributes
- ---------
- tms : fwk.Timers object
- dictionary of timers
- msh : tbox.MshData object
- mesh data structure
- mshWriter : tbox.MshExport object
- mesh data writer
- pbl : dart.Probelm object
- problem formulation
- bnd : array of dart.Body objects
- bodies of interest immersed in fluid
- solver : dart.Solver object
- full potential solver
- """
- def __init__(self, p):
- """Setup and configure
-
- Parameters
- ----------
- p : dict
- dictionary of parameters to configure DART
- """
- # timer
- self.tms = fwk.Timers()
- self.tms["total"].start()
-
- # basic checks
- if p['Dim'] != 2 and p['Dim'] != 3:
- raise RuntimeError('Problem dimension should be 2 or 3, but ' + p['Dim'] + ' was given!\n')
- if p['Dim'] == 2:
- if 'AoS' in p and p['AoS'] != 0:
- raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
- beta = 0
- if 'Symmetry' in p:
- raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
- else:
- if 'AoS' in p:
- if p['AoS'] != 0 and 'Symmetry' in p:
- raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
- beta = p['AoS']*math.pi/180
- else:
- beta = 0
- # basic config
- if p['Format'] == 'vtk':
- try:
- import tboxVtk
- Writer = tboxVtk.VtkExport
- print("Found VTK libraries! 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")
-
- # mesh the airfoil
- print(ccolors.ANSI_BLUE + 'Meshing the airfoil...' + ccolors.ANSI_RESET)
- self.msh = gmsh.MeshLoader(p['File']).execute(**p['Pars'])
- if p['Dim'] == 2:
- mshCrck = tbox.MshCrack(self.msh, p['Dim'])
- mshCrck.setCrack(p['Wake'])
- mshCrck.addBoundaries([p['Fluid'], p['Farfield'][-1], p['Wing']])
- mshCrck.run()
- else:
- for i in range(0, len(p['Wakes'])):
- mshCrck = tbox.MshCrack(self.msh, p['Dim'])
- mshCrck.setCrack(p['Wakes'][i])
- mshCrck.addBoundaries([p['Fluid'], p['Farfield'][-1], p['Wings'][i]])
- if 'Fuselage' in p:
- mshCrck.addBoundaries([p['Fuselage']])
- if 'Symmetry' in p:
- mshCrck.addBoundaries([p['Symmetry']])
- mshCrck.setExcluded(p['WakeTips'][i])
- mshCrck.run()
- tbox.GmshExport(self.msh).save(self.msh.name)
- del mshCrck
- self.mshWriter = Writer(self.msh)
- print('\n')
-
- # initialize the problem
- self.pbl = dart.Problem(self.msh, p['Dim'], 0., beta, p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref'])
- # add a fluid "air"
- self.pbl.set(dart.Fluid(self.msh, p['Fluid'], p['M_inf'], p['Dim'], 0., beta))
- # add initial condition
- self.pbl.set(dart.Initial(self.msh, p['Fluid'], p['Dim'], 0., beta))
- # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
- for i in range (0, len(p['Farfield'])):
- self.pbl.add(dart.Freestream(self.msh, p['Farfield'][i], p['Dim'], 0., beta))
- # add solid boundaries
- if p['Dim'] == 2:
- self.bnd = dart.Body(self.msh, [p['Wing'], p['Fluid']])
- self.pbl.add(self.bnd)
- else:
- self.bnd = []
- for bd in p['Wings']:
- bnd = dart.Body(self.msh, [bd, p['Fluid']])
- self.bnd.append(bnd)
- self.pbl.add(bnd)
- if 'Fuselage' in p:
- bnd = dart.Body(self.msh, [p['Fuselage'], p['Fluid']])
- self.bnd.append(bnd)
- self.pbl.add(bnd)
- # add Wake/Kutta boundary conditions
- if p['Dim'] == 2:
- self.pbl.add(dart.Wake(self.msh, [p['Wake'], p['Wake']+'_', p['Fluid']]))
- self.pbl.add(dart.Kutta(self.msh, [p['Te'], p['Wake']+'_', p['Wing'], p['Fluid']]))
- else:
- for i in range(0, len(p['Wakes'])):
- self.pbl.add(dart.Wake(self.msh, [p['Wakes'][i], p['Wakes'][i]+'_', p['Fluid'], p['TeTips'][i]]))
-
- # initialize the linear (inner) solver
- print(ccolors.ANSI_BLUE + 'Initializing the solver...' + ccolors.ANSI_RESET)
- if p['LSolver'] == 'Pardiso':
- linsol = LinearSolver().pardiso()
- elif p['LSolver'] == 'GMRES':
- linsol = tbox.Gmres()
- linsol.setFillFactor(p['G_fill'])
- linsol.setRestart(p['G_restart'])
- linsol.setTolerance(p['G_tol'])
- elif p['LSolver'] == 'MUMPS':
- linsol = LinearSolver().mumps()
- elif p['LSolver'] == 'SparseLU':
- linsol = tbox.SparseLu()
- 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, 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, 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')
diff --git a/dart/api/internal/polar.py b/dart/api/internal/polar.py
index 75644e5..69f5968 100644
--- a/dart/api/internal/polar.py
+++ b/dart/api/internal/polar.py
@@ -19,19 +19,22 @@
## Polar
# Adrien Crovato
#
-# Compute lift and polar curves
+# Compute aerodynamic load coefficients as a function of angle of attack
import math
import dart.utils as dU
import tbox.utils as tU
+import fwk
from fwk.coloring import ccolors
-from .config import Config
+from ..core import initDart
-class Polar(Config):
+class Polar:
"""Utility to compute the polar of a lifting surface
Attributes
----------
+ tms : fwk.Timers object
+ dictionary of timers
dim : int
problem dimensions
format : str
@@ -59,8 +62,11 @@ class Polar(Config):
p : dict
dictionary of parameters to configure DART
"""
+ # basic init
+ self.tms = fwk.Timers()
+ self.tms["total"].start()
+ self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis')
# save some parameters for later use
- self.dim = p['Dim']
self.format = p['Format']
try:
self.slice = p['Slice']
@@ -70,8 +76,6 @@ class Polar(Config):
self.tag = None
self.mach = p['M_inf']
self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
- # basic init
- Config.__init__(self, p)
def run(self):
"""Compute the polar
@@ -90,20 +94,19 @@ class Polar(Config):
# run the solver and save the results
print(ccolors.ANSI_BLUE + 'Running the solver for', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET)
self.tms["solver"].start()
- self.solver.run()
+ self.sol.run()
self.tms["solver"].stop()
- self.solver.save(self.mshWriter, ac)
- print('\n')
+ self.sol.save(self.wrtr, ac)
# extract Cp
if self.dim == 2:
- Cp = dU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
+ Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
tU.write(Cp, 'Cp_airfoil_'+str(ac)+'.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
elif self.dim == 3 and self.format == 'vtk' and self.slice:
dU.writeSlices(self.msh.name, self.slice, self.tag, ac)
# extract force coefficients
- self.Cls.append(self.solver.Cl)
- self.Cds.append(self.solver.Cd)
- self.Cms.append(self.solver.Cm)
+ self.Cls.append(self.sol.Cl)
+ self.Cds.append(self.sol.Cd)
+ self.Cms.append(self.sol.Cm)
# end loop
ac+=1
diff --git a/dart/api/internal/trim.py b/dart/api/internal/trim.py
index 34b3d0b..2555218 100644
--- a/dart/api/internal/trim.py
+++ b/dart/api/internal/trim.py
@@ -24,15 +24,18 @@
import math
import dart.utils as dU
import tbox.utils as tU
+import fwk
from fwk.coloring import ccolors
-from .config import Config
+from ..core import initDart
-class Trim(Config):
+class Trim:
"""Utility to perform the trim analysis of a given lifting configuration
Find the angle of attack to reach a desired lift coefficient
Attributes
----------
+ tms : fwk.Timers object
+ dictionary of timers
dim : int
problem dimensions
format : str
@@ -51,8 +54,18 @@ class Trim(Config):
initial (estimated) slope of lift curve
"""
def __init__(self, p):
+ """Setup and configure
+
+ Parameters
+ ----------
+ p : dict
+ dictionary of parameters to configure DART
+ """
+ # basic init
+ self.tms = fwk.Timers()
+ self.tms["total"].start()
+ self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis')
# save some parameters for later use
- self.dim = p['Dim']
self.format = p['Format']
try:
self.slice = p['Slice']
@@ -64,8 +77,6 @@ class Trim(Config):
self.clTgt = p['CL']
self.alpha = p['AoA']*math.pi/180
self.dCl = p['dCL']
- # basic init
- Config.__init__(self, p)
def run(self):
"""Perform the trim analysis
@@ -79,26 +90,25 @@ class Trim(Config):
self.pbl.update(self.alpha)
# run the solver
self.tms["solver"].start()
- success = self.solver.run()
+ success = self.sol.run()
if not success:
raise RuntimeError('Flow solver diverged!\n')
self.tms["solver"].stop()
# update slope
if it != 0:
- self.dCl = (self.solver.Cl - Cl_) / (self.alpha - alpha_)
+ self.dCl = (self.sol.Cl - Cl_) / (self.alpha - alpha_)
# break or store values and update AoA
- if abs(self.clTgt - self.solver.Cl) < 0.005 or math.isnan(self.solver.Cl):
+ if abs(self.clTgt - self.sol.Cl) < 0.005 or math.isnan(self.sol.Cl):
break
else:
- Cl_ = self.solver.Cl
+ Cl_ = self.sol.Cl
alpha_ = self.alpha
- dAlpha = (self.clTgt - self.solver.Cl) / self.dCl
+ dAlpha = (self.clTgt - self.sol.Cl) / self.dCl
self.alpha += dAlpha
it += 1
# save results
- self.solver.save(self.mshWriter)
- print('\n')
+ self.sol.save(self.wrtr)
# extract Cp
if self.dim == 2:
Cp = dU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
@@ -107,12 +117,12 @@ class Trim(Config):
dU.writeSlices(self.msh.name, self.slice, self.tag)
def disp(self):
- """DIsplay the results
+ """Display the results
"""
# display results
print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET)
print(" M alpha dCl Cl Cd Cm")
- print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.solver.Cl, self.solver.Cd, self.solver.Cm))
+ print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.sol.Cl, self.sol.Cd, self.sol.Cm))
print('\n')
# display timers
self.tms["total"].stop()
diff --git a/dart/api/mphys.py b/dart/api/mphys.py
index 06f0ca0..2edba1e 100644
--- a/dart/api/mphys.py
+++ b/dart/api/mphys.py
@@ -433,14 +433,10 @@ class DartBuilder(Builder):
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
@@ -453,192 +449,15 @@ class DartBuilder(Builder):
Parameters
----------
params : dict
- Dictonnary of parameters for the mesh, solver, ...
+ Dictonnary of parameters to configure DART
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['Q_inf']
- 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 body of interest (FSI/OPT) is always first given body
- 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.Fluid(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.Body(self.__msh, [params['Wing'], params['Fluid']])
- self.__bnd = bnd
- self.__pbl.add(bnd)
- else:
- self.__bnd = None
- for bd in params['Wings']:
- bnd = dart.Body(self.__msh, [bd, params['Fluid']])
- if self.__bnd is None:
- self.__bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
- self.__pbl.add(bnd)
- if 'Fuselage' in params:
- bnd = dart.Body(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
+ from .core import initDart
+ self.__dim, self.__qinf, _, self.__wrtr, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(params, scenario=scenario, task=task)
print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
def get_mesh_coordinate_subsystem(self):
diff --git a/dart/cases/coyote.py b/dart/cases/coyote.py
index e23edac..6dda2e0 100644
--- a/dart/cases/coyote.py
+++ b/dart/cases/coyote.py
@@ -25,7 +25,12 @@
def getParam():
import os.path
+ from fwk.wutils import parseargs
p = {}
+ # Arguments
+ args = parseargs()
+ p['Threads'] = args.k
+ p['Verb'] = args.verb
# Specific
scl = 2 # scaling factor for lifting surface mesh size
wrtems = scl*0.036 # wing root trailing edge mesh size
@@ -61,17 +66,13 @@ def getParam():
p['y_ref'] = 1.8 # Reference point for moment computation (y)
p['z_ref'] = 0.8 # Reference point for moment computation (z)
# Numerical
- p['LSolver'] = 'GMRES' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
+ p['LSolver'] = 'GMRES' # Linear solver (PARDISO, GMRES, or MUMPS)
p['G_fill'] = 2 # Fill-in factor for GMRES preconditioner
p['G_tol'] = 1e-5 # Tolerance for GMRES
p['G_restart'] = 50 # Restart for GMRES
- p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 15 # Solver maximum number of iterations
- p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
- p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
- p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
return p
def main():
diff --git a/dart/cases/n0012.py b/dart/cases/n0012.py
index f984557..ed9f3da 100644
--- a/dart/cases/n0012.py
+++ b/dart/cases/n0012.py
@@ -21,7 +21,12 @@
def getParam():
import os.path
+ from fwk.wutils import parseargs
p = {}
+ # Arguments
+ args = parseargs()
+ p['Threads'] = args.k
+ p['Verb'] = args.verb
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012.geo') # Input file containing the model
p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model
@@ -45,12 +50,10 @@ def getParam():
p['y_ref'] = 0. # Reference point for moment computation (y)
p['z_ref'] = 0. # Reference point for moment computation (z)
# Numerical
- p['LSolver'] = 'SparseLU' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
- p['NSolver'] = 'Picard' # Noninear solver type (Picard or Newton)
+ p['LSolver'] = 'GMRES' # Linear solver (PARDISO, GMRES, or MUMPS)
p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 20 # Solver maximum number of iterations
- p['Relaxation'] = 0.7 # Relaxation parameter for Picard solver
return p
def main():
diff --git a/dart/cases/n0012_3.py b/dart/cases/n0012_3.py
index 1133860..57700c2 100644
--- a/dart/cases/n0012_3.py
+++ b/dart/cases/n0012_3.py
@@ -21,9 +21,14 @@
def getParam():
import os.path
+ from fwk.wutils import parseargs
p = {}
# Specific
span = 1. # span length
+ # Arguments
+ args = parseargs()
+ p['Threads'] = args.k
+ p['Verb'] = args.verb
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012_3.geo') # Input file containing the model
p['Pars'] = {'spn' : span, 'lgt' : 3., 'wdt' : 3., 'hgt' : 3., 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.} # Parameters for input file model
@@ -51,14 +56,10 @@ def getParam():
p['y_ref'] = 0. # Reference point for moment computation (y)
p['z_ref'] = 0. # Reference point for moment computation (z)
# Numerical
- p['LSolver'] = 'MUMPS' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
- p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
+ p['LSolver'] = 'MUMPS' # Linear solver (PARDISO, GMRES, or MUMPS)
p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 10 # Solver maximum number of iterations
- p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
- p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
- p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
return p
def main():
diff --git a/dart/cases/n64A410.py b/dart/cases/n64A410.py
index 0ef97fc..d033dcb 100644
--- a/dart/cases/n64A410.py
+++ b/dart/cases/n64A410.py
@@ -21,7 +21,12 @@
def getParam():
import os.path
+ from fwk.wutils import parseargs
p = {}
+ # Arguments
+ args = parseargs()
+ p['Threads'] = args.k
+ p['Verb'] = args.verb
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n64A410.geo') # Input file containing the model
p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model
@@ -45,14 +50,10 @@ def getParam():
p['y_ref'] = 0. # Reference point for moment computation (y)
p['z_ref'] = 0. # Reference point for moment computation (z)
# Numerical
- p['LSolver'] = 'Pardiso' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
- p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
+ p['LSolver'] = 'PARDISO' # Linear solver (PARDISO, GMRES, or MUMPS)
p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 25 # Solver maximum number of iterations
- p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
- p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
- p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
return p
def main():
diff --git a/dart/cases/rae2822.py b/dart/cases/rae2822.py
index cf51570..1b9c6cb 100644
--- a/dart/cases/rae2822.py
+++ b/dart/cases/rae2822.py
@@ -21,7 +21,12 @@
def getParam():
import os.path
+ from fwk.wutils import parseargs
p = {}
+ # Arguments
+ args = parseargs()
+ p['Threads'] = args.k
+ p['Verb'] = args.verb
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/rae2822.geo') # Input file containing the model
p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model
@@ -45,17 +50,13 @@ def getParam():
p['y_ref'] = 0. # Reference point for moment computation (y)
p['z_ref'] = 0. # Reference point for moment computation (z)
# Numerical
- p['LSolver'] = 'GMRES' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
+ p['LSolver'] = 'GMRES' # Linear solver (PARDISO, GMRES, or MUMPS)
p['G_fill'] = 2 # Fill-in factor for GMRES preconditioner
p['G_tol'] = 1e-5 # Tolerance for GMRES
p['G_restart'] = 50 # Restart for GMRES
- p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 10 # Solver maximum number of iterations
- p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
- p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
- p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
return p
def main():
diff --git a/dart/cases/wbht.py b/dart/cases/wbht.py
index 4255fff..4744d14 100644
--- a/dart/cases/wbht.py
+++ b/dart/cases/wbht.py
@@ -24,6 +24,7 @@
def getParam():
import os.path
+ from fwk.wutils import parseargs
p = {}
# Specific
scl = 2 # scaling factor for lifting surface mesh size
@@ -40,6 +41,10 @@ def getParam():
fums = 0.15 # fuselage mesh size
fucms = 0.01 # fuselage caps mesh size
fams = 5. # farfield mesh size
+ # Arguments
+ args = parseargs()
+ p['Threads'] = args.k
+ p['Verb'] = args.verb
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/wbht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
@@ -73,17 +78,13 @@ def getParam():
p['y_ref'] = 0.66 # Reference point for moment computation (y)
p['z_ref'] = -0.5 # Reference point for moment computation (z)
# Numerical
- p['LSolver'] = 'GMRES' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
+ p['LSolver'] = 'GMRES' # Linear solver (PARDISO, GMRES, or MUMPS)
p['G_fill'] = 2 # Fill-in factor for GMRES preconditioner
p['G_tol'] = 1e-5 # Tolerance for GMRES
p['G_restart'] = 50 # Restart for GMRES
- p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 15 # Solver maximum number of iterations
- p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
- p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
- p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
return p
def main():
diff --git a/dart/cases/wht.py b/dart/cases/wht.py
index c5f90f5..ce73017 100644
--- a/dart/cases/wht.py
+++ b/dart/cases/wht.py
@@ -21,6 +21,7 @@
def getParam():
import os.path
+ from fwk.wutils import parseargs
p = {}
# Specific
scl = 2 # scaling factor for lifting surface mesh size
@@ -35,6 +36,10 @@ def getParam():
ttlems = scl*0.0058 # tail tip leading mesh size
tttems = scl*0.0058 # tail tip trailing mesh size
fams = 5. # farfield mesh size
+ # Arguments
+ args = parseargs()
+ p['Threads'] = args.k
+ p['Verb'] = args.verb
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/wht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
@@ -67,14 +72,10 @@ def getParam():
p['y_ref'] = 0.0 # Reference point for moment computation (y)
p['z_ref'] = 0.0 # Reference point for moment computation (z)
# Numerical
- p['LSolver'] = 'Pardiso' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
- p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
+ p['LSolver'] = 'PARDISO' # Linear solver (PARDISO, GMRES, or MUMPS)
p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 15 # Solver maximum number of iterations
- p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
- p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
- p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
return p
def main():
--
GitLab