From 6becdc60f493a282580ced9ec718bbecfa99ff9c Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Tue, 21 Jan 2025 14:52:46 +0100
Subject: [PATCH] (api) Fix api for 3D computations

Plus fixed a bug that would set span to 0
---
 blast/api/blaster_api.py                    |  11 +-
 blast/tests/{test_api.py => test_api_2D.py} |  12 +-
 blast/tests/test_api_3D.py                  | 161 ++++++++++++++++++++
 blast/utils.py                              |   7 +-
 4 files changed, 177 insertions(+), 14 deletions(-)
 rename blast/tests/{test_api.py => test_api_2D.py} (95%)
 create mode 100644 blast/tests/test_api_3D.py

diff --git a/blast/api/blaster_api.py b/blast/api/blaster_api.py
index 3d34d6d..5551c31 100644
--- a/blast/api/blaster_api.py
+++ b/blast/api/blaster_api.py
@@ -72,7 +72,9 @@ def init_blaster(cfg, icfg, iSolverName='DART', task='analysis'):
     _verbose = cfg.get('Verb', 0)
     _xtrF = cfg.get('xtrF', [-1, -1])
     _span = cfg.get('span', 1.)
-    _nSections = cfg.get('nSections', 1)
+    if 'sections' not in cfg:
+        raise RuntimeError('BLASTERAPI: Missing sections in the configuration')
+    _nSections = len(cfg['sections'])
 
     if _Minf <= 0:
         raise RuntimeError('BLASTERAPI: Invalid Mach number', _Minf)
@@ -106,6 +108,13 @@ def init_blaster(cfg, icfg, iSolverName='DART', task='analysis'):
     # Coupler
     coupler = Coupler(solverAPI, vSolver, _maxCouplIter = _couplIter, _couplTol = _couplTol, _iterPrint = _iterPrint, _resetInv = _resetInv)
 
+    # Sanity check
+    if solverAPI.getnDim() == 3:
+        if 'span' not in cfg:
+            raise RuntimeError('BLASTERAPI: Missing span for 3D computation')
+        if cfg['interpolator'] != 'Rbf':
+            raise RuntimeError('BLASTERAPI: 3D computation only supports Rbf interpolator')
+
     # Adjoint objects
     cAdj = None
     if task == 'optimization':
diff --git a/blast/tests/test_api.py b/blast/tests/test_api_2D.py
similarity index 95%
rename from blast/tests/test_api.py
rename to blast/tests/test_api_2D.py
index ef0f793..a9530b5 100644
--- a/blast/tests/test_api.py
+++ b/blast/tests/test_api_2D.py
@@ -18,25 +18,15 @@
 
 # @author Paul Dechamps
 # @date 2024
-# Test the blaster adjoint implementation.
+# Test the blaster api for a 2D problem.
 
 # 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 vutils
 
 def cfgInviscid(nthrds, verb):
     import os.path
diff --git a/blast/tests/test_api_3D.py b/blast/tests/test_api_3D.py
new file mode 100644
index 0000000..6959252
--- /dev/null
+++ b/blast/tests/test_api_3D.py
@@ -0,0 +1,161 @@
+#!/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 api for a 2D problem.
+
+# Imports.
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+import blast.utils as vutils
+
+import math
+import numpy as np
+import os
+
+from blast.api.blaster_api import init_blaster
+
+def cfgInviscid(nthrds, verb):
+    # Parameters
+    return {
+    # Options
+    'scenario' : 'aerodynamic',
+    'task' : 'analysis',
+    'Threads' : nthrds, # number of threads
+    'Verb' : verb, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.abspath(os.path.join(os.path.abspath(__file__), '../../models/dart/rae_3.geo')), # Input file containing the model
+    'Pars' : {'spn' : 1.0, 'lgt' : 6.0, 'wdt' : 3.0, 'hgt' : 6.0, 'msN' : 0.02, 'msF' : 1}, # parameters for input file model
+    'Dim' : 3, # problem dimension
+    'Format' : 'vtk', # 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)
+    'Wings' : ['wing'], # LIST of names of physical groups containing the lifting surface boundary
+    'Wakes' : ['wake'], # LIST of names of physical group containing the wake
+    'WakeTips' : ['wakeTip'], # LIST of names of physical group containing the edge of the wake
+    'Tes' : ['te'], # LIST of names of physical group containing the trailing edge
+    'Symmetry': 'symmetry', # name of physical group containing the symmetry BC
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : 0.8, # freestream Mach number
+    'AoA' : 0.0, # freestream angle of attack
+    # Geometry
+    'S_ref' : 1.0, # reference surface length
+    'c_ref' : 1.0, # reference chord length
+    'x_ref' : 0.0, # reference point for moment computation (x)
+    'y_ref' : 0.0, # reference point for moment computation (y)
+    'z_ref' : 0.0, # reference point for moment computation (z)
+    # Numerical
+    'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+    'G_fill' : 2, # fill-in factor for GMRES preconditioner
+    'G_tol' : 1e-5, # tolerance for GMRES
+    'G_restart' : 50, # restart for GMRES
+    'Rel_tol' : 1e-6, # relative tolerance on solver residual
+    'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+    'Max_it' : 50 # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 1e7,       # Freestream Reynolds number
+        'Minf' : 0.8,     # Freestream Mach number (used for the computation of the time step only)
+        'CFL0' : 1,       # Inital CFL number of the calculation
+
+        'sections' : np.linspace(0.01, 0.95, 3),    # Sections on the wing
+        'writeSections': [0.2, 0.4, 0.6, 0.8, 1.0], # Spanwise locations to write the solution
+        'Sym':[0.],                                 # Symmetry plane
+        'span': 1.,                                 # Span of the wing
+        'interpolator': 'Rbf',                      # Interpolator type
+        'rbftype': 'linear',                        # Radial basis function (rbf) type
+        'smoothing': 1e-8,                          # rbf smoothing factor
+        'degree': 0,                                # Degree of the interpolator
+        'neighbors': 10,                            # Number of neighbors for the interpolator
+        'saveTag': 4,                               # Tag to save the solution with VTK
+
+        'Verb': verb,       # Verbosity level of the solver
+        'couplIter': 5,     # Maximum number of iterations
+        'couplTol' : 5e-2,  # Tolerance of the VII methodology
+        'iterPrint': 1,     # int, number of iterations between outputs
+        'resetInv' : False, # bool, flag to reset the inviscid calculation at every iteration.
+        'xtrF' : [0., 0.],  # Forced transition locations
+    }
+
+def main():
+    # Timer.
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+
+    parsViscous = {'spn': icfg['Pars']['spn'], 'lgt': icfg['Pars']['lgt'],
+                   'nLe': 25, 'nMid': 50, 'nTe': 10, 'nSpan': 60, 'nWake': 25,
+                   'progLe': 1.1, 'progMid': 1.0,  'progTe': 1.0, 'progSpan': 1.0, 'progWake': 1.15}
+    
+    vMeshFile = os.path.abspath(os.path.join(os.path.abspath(__file__), '../../models/dart/rae_3_visc.geo'))
+    vMsh = vutils.mesh(vMeshFile, parsViscous)
+    vcfg['vMsh'] = vMsh
+
+    # Init via API
+    obj = init_blaster(vcfg, icfg, iSolverName='DART', task='analysis')
+    coupler, isol, vsol = obj['coupler'], obj['isol'], obj['vsol']
+
+    # Init via utils
+    coupler2, isol2, vsol2 = vutils.initBlast(icfg, vcfg)
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    aeroCoeffs = coupler.run()
+    aeroCoeffs2 = coupler2.run()
+    tms['solver'].stop()
+
+    # Display results.
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd  Cd_wake      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f} {8:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdf + isol.getCd(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+
+    # Save results
+    isol.save(sfx='_viscous')
+    tms['total'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+    
+    # Test solution
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    tests.add(CTest('Cl', isol.getCl(), isol2.getCl(), 1e-12))
+    tests.add(CTest('Cd wake', vsol.Cdt, vsol2.Cdt, 1e-12, forceabs=True))
+    tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), vsol2.Cdf + isol2.getCd(), 1e-12, forceabs=True))
+    tests.add(CTest('Cdf', vsol.Cdf, vsol2.Cdf, 1e-12, forceabs=True))
+    tests.add(CTest('Iterations', len(aeroCoeffs), len(aeroCoeffs2), 0, forceabs=True))
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()
diff --git a/blast/utils.py b/blast/utils.py
index d863f53..8bb1cb5 100644
--- a/blast/utils.py
+++ b/blast/utils.py
@@ -3,7 +3,7 @@ from fwk.wutils import parseargs
 from fwk.coloring import ccolors
 import numpy as np
 
-def initBL(Re, Minf, CFL0, nSections, xtrF = [None, None], span=0, verb=None):
+def initBL(Re, Minf, CFL0, nSections, span, xtrF = [None, None], verb=None):
     """ Initialize boundary layer solver.
     
     Params
@@ -30,6 +30,9 @@ def initBL(Re, Minf, CFL0, nSections, xtrF = [None, None], span=0, verb=None):
     if nSections < 0:
         print(ccolors.ANSI_RED, "blast::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET)
         raise RuntimeError("Invalid parameter")
+    if span <= 0.:
+        print(ccolors.ANSI_RED, "blast::vUtils Fatal error : Negative span.", span, ccolors.ANSI_RESET)
+        raise RuntimeError("Invalid parameter")
     if verb is None:
       args = parseargs()
       verbose = args.verb
@@ -70,7 +73,7 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
     if 'sfx' not in vconfig:
         vconfig['sfx'] = ''
     # Viscous solver
-    vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'], verb=vconfig['Verb'])
+    vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], span=vconfig['span'], xtrF=vconfig['xtrF'], verb=vconfig['Verb'])
 
     # Inviscid solver
     if iSolver == 'DART':
-- 
GitLab