(feat) Update API

API only works for 2D cases

blast/api/core.py → blast/api/blaster_api.py

@@ -19,7 +19,9 @@
 ## Initialize blast computation
 # Paul Dechamps
 
-def initBlast(cfg, icfg, iSolverName='DART'):
+AVAILABLE_SOLVERS = ['DART']
+
+def init_blaster(cfg, icfg, iSolverName='DART', task='analysis'):
     """
     Inputs
     ------
@@ -55,60 +57,61 @@ def initBlast(cfg, icfg, iSolverName='DART'):
     from blast.coupler import Coupler
     import blast
 
-    if iSolverName == 'DART':
+    if iSolverName not in AVAILABLE_SOLVERS:
+        raise RuntimeError('BLASTERAPI: Invalid inviscid solver name', iSolverName)
+    elif iSolverName == 'DART':
         from blast.interfaces.dart.DartInterface import DartInterface as interface
 
-    # Check viscous solver parameters.
-    if 'Re' in cfg and cfg['Re'] > 0:
-        _Re = cfg['Re']
-    else:
-        raise RuntimeError('Missing or invalid Reynolds number')
-    if 'Minf' in cfg and cfg['Minf'] > 0:
-        _Minf = cfg['Minf']
-    else:
-        _Minf = 0.1
-    if 'CFL0' in cfg and cfg['CFL0'] > 0:
-        _CFL0 = cfg['CFL0']
-    else:
-        _CFL0 = 1
-    if 'Verb' in cfg and 0<= cfg['Verb'] <= 3:
-        _verbose = cfg['Verb']
-    else:
-        _verbose = 1
-    _xtrF = [-1, -1]
-    if 'xtrF' in cfg:
-        for i, xtr in enumerate(cfg['xtrF']):
-            if not(0 <= xtr <= 1) and xtr != -1:
-                raise RuntimeError("Incorrect forced transition location.")
-            _xtrF[i] = xtr
-    _span = 0
-    _nSections = 1 
-
-    # Check coupler parameters.
-    if 'couplIter' in cfg and cfg['couplIter'] > 0:
-        __couplIter = cfg['couplIter']
-    else:
-        __couplIter = 150
-    if 'couplTol' in cfg:
-        __couplTol = cfg['couplTol']
-    else:
-        __couplTol = 1e-4
-    if 'iterPrint' in cfg:
-        __iterPrint = cfg['iterPrint']
-    else:
-        __iterPrint = 1
-    if 'resetInv' in cfg:
-        __resetInv = cfg['resetInv']
-    else:
-        __resetInv = False
+    # Viscous solver
+    if 'Re' not in cfg or cfg['Re'] <= 0.:
+        raise RuntimeError('BLASTERAPI: Missing or invalid Reynolds number')
+    _Re = cfg['Re']
+
+    _Minf = cfg.get('Minf', 0.1)
+    _cfl0 = cfg.get('CFL0', 1.)
+    _verbose = cfg.get('Verb', 0)
+    _xtrF = cfg.get('xtrF', [-1, -1])
+    _span = cfg.get('span', 1.)
+    _nSections = cfg.get('nSections', 1)
+
+    if _Minf <= 0:
+        raise RuntimeError('BLASTERAPI: Invalid Mach number', _Minf)
+    if _cfl0 <= 0:
+        raise RuntimeError('BLASTERAPI: Invalid CFL number', _cfl0)
+    for xtr in _xtrF:
+        if not(0 <= xtr <= 1) and xtr != -1:
+            raise RuntimeError("BLASTERAPI: Incorrect forced transition location.", xtr)
+    if _span < 0:
+        raise RuntimeError('BLASTERAPI: Invalid span', _span)
+    if _nSections <= 0:
+        raise RuntimeError('BLASTERAPI: Invalid number of sections', _nSections)
+
+    # Coupler
+    _couplIter = cfg.get('couplIter', 150)
+    _couplTol = cfg.get('couplTol', 1e-4)
+    _iterPrint = cfg.get('iterPrint', 1)
+    _resetInv = cfg.get('resetInv', False)
+
+    if _couplIter < 0:
+        raise RuntimeError('BLASTERAPI: Invalid number of coupling iterations', _couplIter)
+    if _couplTol <= 0:
+        raise RuntimeError('BLASTERAPI: Invalid coupling tolerance', _couplTol)
+    if _iterPrint < 0:
+        raise RuntimeError('BLASTERAPI: Invalid iteration print frequency', _iterPrint)
 
     # Viscous solver object.
-    vSolver = blast.Driver(_Re, _Minf, _CFL0, _nSections, _xtrF[0], _xtrF[1], _span, _verbose =_verbose)
+    vSolver = blast.Driver(_Re, _Minf, _cfl0, _nSections, _xtrF[0], _xtrF[1], _span, _verbose =_verbose)
     # Solvers interface.
-    solversAPI = interface(icfg, vSolver, cfg)
+    solverAPI = interface(icfg, vSolver, cfg)
     # Coupler
-    coupler = Coupler(solversAPI, vSolver, _maxCouplIter = __couplIter, _couplTol = __couplTol, _iterPrint = __iterPrint, _resetInv = __resetInv)
+    coupler = Coupler(solverAPI, vSolver, _maxCouplIter = _couplIter, _couplTol = _couplTol, _iterPrint = _iterPrint, _resetInv = _resetInv)
+
+    # Adjoint objects
+    cAdj = None
+    if task == 'optimization':
+        cAdj = blast.CoupledAdjoint(solverAPI.adjointSolver, vSolver)
 
     return {'coupler': coupler,
-            'solversAPI': solversAPI,
-            'vSolver': vSolver}
+            'isol': solverAPI,
+            'vsol': vSolver,
+            'adj': cAdj}

blast/tests/apiTest.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2022 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+# @author Paul Dechamps
+# @date 2024
+# Test the blaster adjoint implementation.
+
+# Imports.
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+import numpy as np
+import os
+
+from matplotlib import pyplot as plt
+from blast.api.blaster_api import init_blaster
+import blast.utils as viscUtils
+
+def cfgInviscid(nthrds, verb):
+    import os.path
+    # Parameters
+    return {
+    # Options
+    'scenario' : 'aerodynamic',
+    'task' : 'optimization',
+    'Threads' : nthrds, # number of threads
+    'Verb' : verb, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.dirname(os.path.abspath(__file__)) + '/../models/dart/n0012.geo', # Input file containing the model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 10, 'msTe' : 0.01, 'msLe' : 0.01}, # parameters for input file model
+    'Dim' : 2, # problem dimension
+    'Format' : 'gmsh', # save format (vtk or gmsh)
+    # Markers
+    'Fluid' : 'field', # name of physical group containing the fluid
+    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+    'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+    'Wake' : 'wake', # LIST of names of physical group containing the wake
+    'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+    'Te' : 'te', # LIST of names of physical group containing the trailing edge
+    'dbc' : False,
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : 0.2, # freestream Mach number
+    'AoA' : 2., # freestream angle of attack
+    # Geometry
+    'S_ref' : 1., # reference surface length
+    'c_ref' : 1., # reference chord length
+    'x_ref' : .25, # reference point for moment computation (x)
+    'y_ref' : 0.0, # reference point for moment computation (y)
+    'z_ref' : 0.0, # reference point for moment computation (z)
+    # Numerical
+    'LSolver' : 'SparseLu', # inner solver (Pardiso, MUMPS or GMRES)
+    'G_fill' : 2, # fill-in factor for GMRES preconditioner
+    'G_tol' : 1e-5, # tolerance for GMRES
+    'G_restart' : 50, # restart for GMRES
+    'Rel_tol' : 1e-6, # relative tolerance on solver residual
+    'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+    'Max_it' : 20, # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 1e6,                     # Freestream Reynolds number
+        'Minf' : 0.2,                   # Freestream Mach number (used for the computation of the time step only)
+        'CFL0' : 1,                     # Inital CFL number of the calculation
+        'Verb': verb,                   # Verbosity level of the solver
+        'couplIter': 100,               # Maximum number of iterations
+        'couplTol' : 1e-2,              # Tolerance of the VII methodology
+        'iterPrint': 20,                # int, number of iterations between outputs
+        'resetInv' : True,              # bool, flag to reset the inviscid calculation at every iteration.
+        'sections' : [0],               # List of sections for boundary layer calculation
+        'xtrF' : [-1, -1],            # Forced transition location
+        'interpolator' : 'Matching',    # Interpolator for the coupling
+    }
+
+def main():
+
+    tms = fwk.Timers()
+
+    tms['total'].start()
+
+    # Parse agrs
+    args = parseargs()
+    verb = args.verb
+    nthrds = args.k
+
+    icfg = cfgInviscid(nthrds, verb)
+    cfg = cfgBlast(verb)
+
+    obj = init_blaster(cfg, icfg, iSolverName='DART', task='optimization')
+    coupler, isol, vsol, adj = obj['coupler'], obj['isol'], obj['vsol'], obj['adj']
+
+    # Forward problem
+    coupler.run()
+
+    # Adjoint problem
+    isol.adjointSolver.run()
+    adj.run()
+
+    vSolution = viscUtils.getSolution(isol.sec, write=False)[0]
+    #plotSolution(vSolution)
+    tms['total'].stop()
+
+    # Make the api crash
+    params = {
+        'Re': -1,
+        'Minf': -1,
+        'CFL0': -1,
+        'couplIter': -1,
+        'couplTol': -1,
+        'iterPrint': -1
+    }
+
+    # Loop over the parameters
+    for param, value in params.items():
+        original = cfg[param]
+        try:
+            cfg[param] = value
+            obj = init_blaster(cfg, icfg, iSolverName='DART', task='optimization')
+            raise AssertionError(f'API initialized with {param} = ', cfg[param])
+        except AssertionError as e:
+            raise RuntimeError(e)
+        except:
+            print(ccolors.ANSI_GREEN + f'API crashed with {param} = ' + str(cfg[param]) + ccolors.ANSI_RESET)
+            cfg[param] = original
+
+    print('Statistics')
+    print(tms)
+    print(ccolors.ANSI_BLUE + 'PyTesting' + ccolors.ANSI_RESET)
+    tests = CTests()
+    tests.add(CTest('Cl', isol.getCl(), 0.226, 1e-2))
+    tests.add(CTest('Cd', isol.getCd()+vsol.Cdf, 0.00587, 1e-3))
+    tests.add(CTest('dCl_dAoA', adj.tdCl_AoA, 5.47, 1e-4))
+    tests.add(CTest('dCd_dAoA', adj.tdCd_AoA, 0.09531, 1e-3))
+    tests.run()
+
+    # eof
+    print('')
+
+def plotSolution(vSolution):
+    plt.figure()
+    plt.plot(vSolution['x'], vSolution['cf'], label='$c_f$')
+    plt.xlim([0, 1])
+    plt.xlabel('x/c')
+    plt.ylabel('$c_f$')
+    plt.legend()
+    plt.grid()
+    plt.savefig('cf.png')
+    plt.draw()
+
+    plt.figure()
+    plt.plot(vSolution['x'], vSolution['ctEq'], label='$ctEq$')
+    plt.xlim([0, 1])
+    plt.xlabel('x/c')
+    plt.ylabel('$ctEq$')
+    plt.legend()
+    plt.grid()
+    plt.savefig('cf.png')
+    plt.draw()
+
+    plt.figure()
+    plt.plot(vSolution['x'], vSolution['theta'], label='$\theta$')
+    plt.plot(vSolution['x'], vSolution['deltaStar'], label='$\delta^*$')
+    plt.plot(vSolution['x'], vSolution['theta']*vSolution['H'], label='$\delta^*from$')
+    plt.xlim([0, 1])
+    plt.xlabel('x/c')
+    plt.ylabel('$\delta^*$')
+    plt.legend()
+    plt.grid()
+    plt.savefig('cf.png')
+    plt.draw()
+
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file