Skip to content
Snippets Groups Projects
  • Paul Dechamps's avatar
    2e2e93bd
    (feat) Added vtk support · 2e2e93bd
    Paul Dechamps authored
    Added vtk support in Dart api. Modified writing of boundary layer slices and introduced new argument of blast - writeSlices - which indicates vtk where to extract cps. The feature is not yet available for the boundary layer slices as this would require interpolation in the viscous solver.
    (feat) Added vtk support
    Paul Dechamps authored
    Added vtk support in Dart api. Modified writing of boundary layer slices and introduced new argument of blast - writeSlices - which indicates vtk where to extract cps. The feature is not yet available for the boundary layer slices as this would require interpolation in the viscous solver.
utils.py 7.36 KiB
import blast
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):
    """ Initialize boundary layer solver.
    
    Params
    ------
    - Re: Flow Reynolds number.
    - Minf: Freestream Mach number.
    - CFL0: Initial CFL number for time integration.
    - nSections: Number of sections in the domain.
    - span: Wing span (not used for 2D computations.
    - verb: Verbosity level of the solver.
    
    Return
    ------
    - solver: blast::Driver class.
    """
    if Re<=0.:
        print(ccolors.ANSI_RED, "blast::vUtils Error : Negative Reynolds number.", Re, ccolors.ANSI_RESET)
        raise RuntimeError("Invalid parameter")
    if Minf<0.:
        print(ccolors.ANSI_RED, "blast::vUtils Error : Negative Mach number.", Minf, ccolors.ANSI_RESET)
        raise RuntimeError("Invalid parameter")
    elif Minf>=1.:
        print(ccolors.ANSI_YELLOW, "blast::vUtils Warning : (Super)sonic freestream Mach number.", Minf, ccolors.ANSI_RESET)
    if nSections < 0:
        print(ccolors.ANSI_RED, "blast::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET)
        raise RuntimeError("Invalid parameter")
    if verb is None:
      args = parseargs()
      verbose = args.verb
    else:
      if not(0<=verb<=3):
        print(ccolors.ANSI_RED, "blast::vUtils Fatal error : verbose not in valid range.", verbose, ccolors.ANSI_RESET)
        raise RuntimeError("Invalid parameter")
      else:
        verbose = verb
    
    for i in range(len(xtrF)):
        if xtrF[i] is None:
            xtrF[i] = -1
        if xtrF[i] != -1 and not(0<= xtrF[i] <= 1):
            raise RuntimeError('Incorrect forced transition location') 

    solver = blast.Driver(Re, Minf, CFL0, nSections, xtrF_top=xtrF[0], xtrF_bot=xtrF[1], _span=span, _verbose=verbose)
    return solver

def initBlast(iconfig, vconfig, iSolver='DART'):
    """Initialize blast coupling objects.
    
    Params
    ------
    - iconfig (dict): Dictionnary to initialize solver 'iSolver'.
    - vconfig (dict): Dictionnary to initialize boundary layer solver.
    - iSolver (string): Name of the inviscid solver to use.

    Return
    ------
    - coupler: coupler object
    - iSolverAPI: api interface of the inviscid solver selected with 'iSolver'.
    - vSolver: blast::Driver class.
    """

    if 'nSections' not in vconfig:
        vconfig['nSections'] = len(vconfig['sections'])
    # Viscous solver
    vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'])

    # Inviscid solver
    if iSolver == 'DART':
        from blast.interfaces.dart.DartInterface import DartInterface as interface
    else:
        print(ccolors.ANSI_RED + 'Solver '+iSolver+' currently not implemented' + ccolors.ANSI_RESET)
        raise RuntimeError
    iSolverAPI = interface(iconfig, vSolver, vconfig)

    # Coupler
    import blast.coupler as blastCoupler
    coupler = blastCoupler.Coupler(iSolverAPI, vSolver, _maxCouplIter=vconfig['couplIter'], _couplTol=vconfig['couplTol'], _iterPrint=vconfig['iterPrint'], _resetInv=vconfig['resetInv'])
    return coupler, iSolverAPI, vSolver

def mesh(file, pars):
    """Initialize mesh and mesh writer

    Parameters
    ----------
    file : str
        file contaning mesh (w/o extention)
    pars : dict
        parameters for mesh
    """
    import tbox.gmsh as tmsh
    # Load the mesh.
    msh = tmsh.MeshLoader(file,__file__).execute(**pars)
    return msh


def getSolution(vSolver, iSec=0):
    """Extract viscous solution.
    Args
    ----
    - vSolver: blast::Driver class. Driver of the viscous calculations
    - iSec (int): Marker of the section (default: 0)
    """
    if iSec<0:
        raise RuntimeError("blast::viscU Invalid section number", iSec)
    
    solverOutput = vSolver.getSolution(iSec)

    sln = { 'theta'    : solverOutput[0],
            'H'        : solverOutput[1],
            'N'        : solverOutput[2],
            'ue'       : solverOutput[3],
            'Ct'       : solverOutput[4],
            'deltaStar': solverOutput[5],
            'Ret'      : solverOutput[6],
            'Cd'       : solverOutput[7],
            'Cf'       : np.asarray(solverOutput[8])*np.asarray(solverOutput[3])**2,
            'Hstar'    : solverOutput[9],
            'Hstar2'   : solverOutput[10],
            'Hk'       : solverOutput[11],
            'Cteq'     : solverOutput[12],
            'Us'       : solverOutput[13],
            'delta'    : solverOutput[14],
            'x'        : solverOutput[15],
            'xoc'      : solverOutput[16],
            'y'        : solverOutput[17],
            'z'        : solverOutput[18],
            'ueInv'    : solverOutput[19],
            'xtrT'     : vSolver.getxtr(iSec, 0),
            'xtrB'     : vSolver.getxtr(iSec, 1),
            'Cdt_w'    : vSolver.Cdt_sec[iSec],
            'Cdf'      : vSolver.Cdf_sec[iSec],
            'Cdp'      : vSolver.Cdp_sec[iSec]
            }
    return sln

def write(wData, Re, toW=['deltaStar', 'H', 'Hstar', 'Cf', 'Ct', 'ue', 'ueInv', 'delta'], sfx=''):
    """Write the results in bl files
    """
    import os
    if not os.path.exists('blSlices'):
        os.makedirs('blSlices')
    # Write
    print('Writing file: /blSlices/bl'+sfx+'.dat...', end = '')
    f = open('blSlices/bl'+sfx+'.dat', 'w+')

    f.write('$Sectional aerodynamic coefficients\n')
    f.write('             Re             Cdw             Cdp             Cdf         xtr_top         xtr_bot\n')
    f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cdt_w'], wData['Cdp'],
                                                                                   wData['Cdf'], wData['xtrT'],
                                                                                   wData['xtrB']))
    f.write('$Boundary layer variables\n')
    f.write('{0:>15s} {1:>15s} {2:>15s} {3:>15s}'.format('x','y', 'z', 'xoc'))


    for s in toW:
        f.write(' {0:>15s}'.format(s))
    f.write('\n')

    for i in range(len(wData['x'])):
        f.write('{0:>15.6f} {1:>15.6f} {2:>15.6f} {3:>15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], (wData['x'][i] - min(wData['x'])) / (max(wData['x']) - min(wData['x']))))
        for s in toW:
            f.write(' {0:15.6f}'.format(wData[s][i]))
        f.write('\n')

    f.close()
    print('done.')

def read(filename, delim=',', skip=1):
    """Read from file and store in data array
    """
    import io
    import numpy as np
    # read file
    fl = io.open(filename, 'r')
    label = fl.readline().split(',')
    fl.close()
    data = np.loadtxt(filename, delimiter=delim, skiprows=skip)
    return data

def plot(cfg):
    from matplotlib import pyplot as plt
    for i in range(len(cfg['curves'])):
        plt.plot(cfg['curves'][i][:,0], cfg['curves'][i][:,1], cfg['ls'][i], color=cfg['color'][i], lw=cfg['lw'][i], label=cfg['labels'][i])
    if 'yreverse' in cfg and cfg['yreverse'] is True: plt.gca().invert_yaxis()
    if 'xreverse' in cfg and cfg['xreverse'] is True: plt.gca().invert_xaxis()
    if 'title' in cfg: plt.title(cfg['title'])
    if 'xlim' in cfg: plt.xlim(cfg['xlim'])
    if 'ylim' in cfg: plt.xlim(cfg['ylim'])
    if 'legend' in cfg and cfg['legend'] is True: plt.legend(frameon=False)
    if 'xlabel' in cfg: plt.xlabel(cfg['xlabel'])
    if 'ylabel' in cfg: plt.ylabel(cfg['ylabel'])
    plt.show()