Implemented HSPM interface

.gitmodules

@@ -4,3 +4,6 @@
 [submodule "modules/UPM"]
 	path = modules/UPM
 	url = git@gitlab.uliege.be:Corentin.Thomee/UPM.git
+[submodule "modules/hspm"]
+	path = modules/hspm
+	url = git@gitlab.uliege.be:Corentin.Thomee/hspm.git
\ No newline at end of file

blast/interfaces/hspm/HSPMInterface.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2024 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.
+
+
+# HSPM interface
+# Corentin Thomee
+
+import numpy as np
+
+from fwk.coloring import ccolors
+from blast.interfaces.DSolversInterface import SolversInterface
+
+class HSPMInterface(SolversInterface):
+    def __init__(self, HSPMCfg, vSolver, _cfg):
+        """
+        Initialize the HSPM solver and interface
+        """
+        try:
+            from modules.hspm.hspm.api.core import initHSPM
+            self.solver = initHSPM(HSPMCfg)
+
+        except:
+            print(ccolors.ANSI_RED, 'Could not load HSPM. Make sure it is installed.', ccolors.ANSI_RESET)
+            raise RuntimeError('Import error')
+
+        SolversInterface.__init__(self, vSolver, _cfg)
+    
+    def run(self):
+        self.solver.solve()
+        return 1 # TODO: implémenter des retours dans la fonction solve en fonction du succès (echec solver linéaire, ITMAX, ...) !!
+    
+    def writeCp(self, sfx=''):
+        self.solver.saveCp('Cp'+sfx+'.dat') # Directly use the HSPM function
+
+    def save(self, sfx='inviscid'):
+        # TODO
+        pass
+
+    def getAoA(self):
+        return self.solver.AoA
+    
+    def getMinf(self):
+        return 0 # Incompressible solver
+    
+    def getCl(self):
+        return self.solver.cl
+
+    def getCd(self):
+        return self.solver.cd
+
+    def getCm(self):
+        return self.solver.cm
+    
+    def getVerbose(self):
+        return 0 # Not necessary
+    
+    def reset(self):
+        pass # Already present in 'run' function in hspm/solver.cpp
+    
+    def imposeInvBC(self, couplIter):
+        """ Extract the inviscid parameters at the edge of the boundary layer
+        and use it as a boundary condition in the viscous solver.
+
+        Params
+        ------
+        - couplIter (int): Coupling iteration.
+        """
+        # Retreive parameters.
+        l = len(self.iBnd[0].nodesCoord[:,3])-1
+        for iRow, row in enumerate(self.iBnd[0].nodesCoord[:,3]):
+            row=int(row)
+
+            # Interpolate from elements to points
+            self.iBnd[0].V[iRow,0] = (self.solver.U[row%l] + self.solver.U[(row-1+l)%l]) / 2
+            self.iBnd[0].V[iRow,1] = (self.solver.V[row%l] + self.solver.V[(row-1+l)%l]) / 2
+            self.iBnd[0].V[iRow,2] = 0 # No z component of speed in 2D
+
+            # Incompressible solver -> M = 0, Rho = cst
+            self.iBnd[0].M[iRow] = 0
+            self.iBnd[0].Rho[iRow] = 1.0
+
+        self.setViscousSolver(couplIter)
+
+    def imposeBlwVel(self):
+        """ Extract the solution of the viscous calculation (blowing velocity)
+        and use it as a boundary condition in the inviscid solver.
+        """
+        self.getBlowingBoundary()
+        # Impose blowing velocities.
+        self.iBnd[0].connectBlowingVel()
+
+        for i, blVel in enumerate(self.iBnd[0].blowingVel):
+            self.solver.imposeBlowingVelocity(i, blVel) # blVel is given on the elements
+
+        # Set the dStar for the inviscid solver
+        """dStarTop = self.vSolver.extractDeltaStar(0, 0)
+        dStarBottom = self.vSolver.extractDeltaStar(0, 1)
+
+        for i in range(len(dStarTop) - 1): # Top surface
+            dStar = (dStarTop[-i-1] + dStarTop[-i-2]) / 2
+            self.solver.setdStar(i, dStar)
+        for i in range(len(dStarBottom) - 1): # Bottom surface
+            dStar = (dStarBottom[i] + dStarBottom[i+1]) / 2
+            self.solver.setdStar(i+len(dStarTop)-1, dStar)"""
+
+    def setMesh(self):
+        if self.cfg['nDim'] == 2:
+            self.getAirfoil()
+        elif self.cfg['nDim'] == 3:
+            raise RuntimeError('HSPMInterface::3D not implemented in UPM')
+        else:
+            raise RuntimeError('HSPMInterface::Could not resolve how to set\
+            the mesh. Either ''nDim'' is incorrect or ''msh'' was not set for\
+            3D computations')
+
+    def getAirfoil(self):
+        """Create data structure for information transfer.
+        """
+        # Number of elements/nodes
+        N_nodes = self.solver.N+1
+        N_elems = self.solver.N
+
+        self.iBnd[0].initStructures(N_nodes) # Initialize the structure for the airfoil boundary
+        
+        connectListNodes = np.array([i for i in range(N_nodes)], dtype=int) # Connectivity list for the nodes
+        connectListElems = np.array([i for i in range(N_elems)], dtype=int) # Connectivity list for the elements
+
+        data = np.zeros((N_nodes, 4))
+        for i in range(N_nodes):
+            data[i,0] = self.solver.x[i]
+            data[i,1] = self.solver.y[i]
+            data[i,2] = 0 # No z component in 2D
+            data[i,3] = i
+        self.iBnd[0].nodesCoord = np.column_stack((data[:,0], data[:,1],\
+                                                   data[:,2], data[:,3]))
+        self.iBnd[0].setConnectList(connectListNodes, connectListElems)

blast/tests/hspm/naca0012_2D.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.
+
+# Imports.
+
+import blast.utils as viscUtils
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+    import os.path
+    # Parameters
+    return {
+	'chord': 1,
+	'cl_tolerance': 1e-5,
+	'it_solver_tolerance': 1e-4,
+	'aoa': 2,
+	'naca':"0012",
+    'N': 100,
+	'filename': os.path.dirname(os.path.abspath(__file__)) + '/../../models/upm/n0012.dat'
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 1e7,       # Freestream Reynolds number
+        'Minf' : 0,     # 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': 500,    # Maximum number of iterations
+        'couplTol' : 1e-4,  # Tolerance of the VII methodology
+        'iterPrint': 1,     # int, number of iterations between outputs
+        'resetInv' : True,  # bool, flag to reset the inviscid calculation at every iteration.
+        'sections' : [0],
+        'hasWake' : False,   # bool, flag to specify if the wake has to be taken into account
+        'xtrF' : [None, None],# Forced transition location
+        'interpolator' : 'Matching',
+        'nDim' : 2
+    }
+
+def main():
+    # Timer.
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+
+    tms['pre'].start()
+    coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg, iSolver='HSPM')
+    tms['pre'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    aeroCoeffs = coupler.run()
+    tms['solver'].stop()
+
+    # Display results.
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, iSolverAPI.getMinf(), iSolverAPI.getAoA()*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
+
+    # Write results to file.
+    vSolution = viscUtils.getSolution(vSolver)
+    vSolution['Cdt_int'] = vSolution['Cdf'] + iSolverAPI.getCd()
+    tms['total'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+
+    # Show results
+    if not args.nogui:
+        iCp = viscUtils.read('Cp_inviscid.dat', delim=" ", skip = 0)
+        vCp = viscUtils.read('Cp_viscous.dat', delim=" ", skip = 0)
+        xfoilCp = viscUtils.read('../../blast/models/upm/references/cp_0012.dat', delim=None, skip = 1)
+        xfoilCpInv = viscUtils.read('../../blast/models/upm/references/cp_0012_inviscid.dat', delim=None, skip = 1)
+        plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
+                            np.column_stack((xfoilCp[:,0], xfoilCp[:,1])),
+                            np.column_stack((iCp[:,0], iCp[:,3])),
+                            np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
+                    'labels': ['Blast (VII)', 'XFOIL VII', 'UPM (inviscid)', 'XFOIL (inviscid)'],
+                    'lw': [3, 3, 2, 2],
+                    'color': ['firebrick', 'darkblue', 'firebrick', 'darkblue'],
+                    'ls': ['-', '-', '--', '--'],
+                    'reverse': True,
+                    'xlim':[0, 1],
+                    'yreverse': True,
+                    'legend': True,
+                    'xlabel': '$x/c$',
+                    'ylabel': '$c_p$'
+                    }
+        viscUtils.plot(plotcp)
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()

blast/utils.py

@@ -76,6 +76,8 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
         from blast.interfaces.dart.DartInterface import DartInterface as interface
     elif iSolver == 'UPM':
         from blast.interfaces.UPM.UPMInterface import UPMInterface as interface
+    elif iSolver == 'HSPM':
+        from blast.interfaces.hspm.HSPMInterface import HSPMInterface as interface
     else:
         print(ccolors.ANSI_RED + 'Solver '+iSolver+' currently not implemented' + ccolors.ANSI_RESET)
         raise RuntimeError

modules/CMakeLists.txt

 ADD_SUBDIRECTORY( dartflo )
-ADD_SUBDIRECTORY( UPM )
\ No newline at end of file
+ADD_SUBDIRECTORY( UPM )
+ADD_SUBDIRECTORY( hspm )
\ No newline at end of file

run.py

@@ -25,12 +25,15 @@ def main():
     fwkdir = os.path.abspath(os.path.join(thisdir, 'modules', 'dartflo', 'ext', 'amfe'))
     dartflodir = os.path.abspath(os.path.join(thisdir, 'modules', 'dartflo'))
     UPMdir = os.path.abspath(os.path.join(thisdir, 'modules', 'UPM'))
+    HSPMdir = os.path.abspath(os.path.join(thisdir, 'modules', 'hspm'))
     if not os.path.isdir(fwkdir):
         raise Exception('blaster/modules/dartflo/amfe not found!\n')
     if not os.path.isdir(dartflodir):
         raise Exception('blaster/modules/dartflo not found!\n')
     if not os.path.isdir(UPMdir):
         raise Exception('blaster/modules/UPM not found!\n')
+    if not os.path.isdir(HSPMdir):
+        raise Exception('blaster/modules/hspm not found!\n')
     sys.path.append(fwkdir)
     sys.path.append(dartflodir)
     sys.path.append(UPMdir)