raeValidation.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 2022

# 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 os.path

import math

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.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/rae_2.geo', # Input file containing the model
    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 10, 'msTe' : 0.0075, 'msLe' : 0.001}, # 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 boundaryk
    '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' : True,
    'Upstream' : 'upstream',
    # Freestream
    'M_inf' : 0.729, # freestream Mach number
    'AoA' : 2.31, # 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' : 50 # solver maximum number of iterations
    }

def cfgBlast(verb):
    return {
        'Re' : 6.5e6,       # Freestream Reynolds number
        'Minf' : 0.73,      # 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': 50,    # Maximum number of iterations
        'couplTol' : 1e-4,  # Tolerance of the VII methodology
        'iterPrint': 5,     # int, number of iterations between outputs
        'resetInv' : True,  # bool, flag to reset the inviscid calculation at every iteration.
        'sections' : [0],
        'xtrF' : [0., 0.],# 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)

    gr = 1.1
    if gr == 1.0:
        icfg['Pars']['msF'] = icfg['Pars']['msLe']
    else:
        n = np.log10(1-(1-gr)*icfg['Pars']['xLgt']/icfg['Pars']['msLe'])/np.log10(gr)
        icfg['Pars']['msF'] = icfg['Pars']['msLe']*gr**(n-1)

    tms['pre'].start()
    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
    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, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))

     # Write results to file.
    vSolution = viscUtils.getSolution(isol.sec)
    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(), 0.759, 5e-2))
    tests.add(CTest('Cd wake', vsol.Cdt, 0.0093, 1e-3, forceabs=True))
    tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0134, 1e-3, forceabs=True))
    tests.add(CTest('Cdf', vsol.Cdf, 0.0068, 1e-3, forceabs=True))
    if icfg['LSolver'] == 'SparseLu':
        tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 29, 0, forceabs=True))
    tests.run()

    expResults = np.loadtxt(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/references/rae2822_AR138_case6.dat')
    expResults[:,1] *= -1

    # Show results
    if not args.nogui:
        iCp = viscUtils.read('Cp_inviscid.dat')
        vCp = viscUtils.read('Cp_viscous.dat')
        plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
                            np.column_stack((iCp[:,0], iCp[:,3])),
                            expResults],
                    'labels': ['Blast (VII)', 'DART (inviscid)', 'Experimental'],
                    'lw': [3, 3, 2, 2],
                    'color': ['darkblue', 'darkblue', 'black'],
                    'ls': ['-', '--', 'o'],
                    'reverse': True,
                    'xlim':[0, 1],
                    'yreverse': True,
                    'legend': True,
                    'xlabel': '$x/c$',
                    'ylabel': '$c_p$'
                    }
        viscUtils.plot(plotcp)

        plotcf = {'curves': [np.column_stack((vSolution[0]['x'], vSolution[0]['cf']))],
                'labels': ['Blast'],
                'lw': [3, 3],
                'color': ['darkblue'],
                'ls': ['-'],
                'reverse': True,
                'xlim':[0, 1],
                'ylim':[0, 0.008],
                'legend': True,
                'xlabel': '$x/c$',
                'ylabel': '$c_f$'
                    }
        viscUtils.plot(plotcf)

    # eof
    print('')

if __name__ == "__main__":
    main()