onera / Fast

Fast CFD solver
https://onera.github.io/Fast/
GNU Lesser General Public License v3.0
14 stars 6 forks source link

GlobalConvergenceHistory is not updated #14

Open Luispain opened 3 months ago

Luispain commented 3 months ago

Hi all,

This message is to inform that GlobalConvergenceHistory does not seem to be updated.

If I run the example script convergenceHistory:

    # - convergenceHistory (pyTree) -
    import Converter.PyTree as C
    import Generator.PyTree as G
    import FastS.PyTree as FastS
    import FastC.PyTree as FastC
    import Connector.PyTree as X
    import Converter.Internal as Internal
    import Initiator.PyTree as I

    ni = 155; dx = 100./(ni-1); dz = 0.01
    a1 = G.cart((-50,-50,0.), (dx,dx,dz), (ni,ni,2))
    a1 = C.fillEmptyBCWith(a1, 'far', 'BCFarfield', dim=2)
    a1 = I.initConst(a1, MInf=0.4, loc='centers')
    a1 = C.addState(a1, 'GoverningEquations', 'Euler')
    a1 = C.addState(a1, MInf=0.4)
    t = C.newPyTree(['Base', a1])

    # Numerics
    modulo_verif = 20
    numb = {}
    numb["temporal_scheme"]    = "implicit"
    numb["ss_iteration"]       = 3
    numb["modulo_verif"]       = modulo_verif
    numz = {}
    numz["time_step"]          = 0.0007
    numz["time_step_nature"]   = "local"
    numz["cfl"]                = 4.0
    numz["scheme"]             = "roe"
    numz["slope"]              = "minmod"
    FastC._setNum2Zones(t, numz) ; FastC._setNum2Base(t, numb)

    (t, tc, metrics) = FastS.warmup(t, None)

    # Number or records to store residuals 
    nrec = 100//modulo_verif

    #To remove old ConvergenceHistory nodes 
    t = C.rmNodes(t, "ZoneConvergenceHistory")

    #Convergence history with nrec records
    FastS.createConvergenceHistory(t, nrec)

    nit = 100; time = 0
    time_step = Internal.getNodeFromName(t, 'time_step')
    time_step = Internal.getValue(time_step)
    for it in range(nit):
        FastS._compute(t, metrics, it, tc)
        if it%modulo_verif == 0:
            FastS.display_temporal_criteria(t, metrics, it)
        time += time_step

    # time stamp
    Internal.createUniqueChild(t, 'Iteration', 'DataArray_t', value=nit)
    Internal.createUniqueChild(t, 'Time', 'DataArray_t', value=time)
    C.convertPyTree2File(t, 'out.cgns')

I can see that GlobalConvergenceHistory node is empty:

image

best regards !

Luis

StephaniePERON commented 2 months ago

Hi Luis, I guess an argument is missing (format='store') in order to store the residuals in the GlobalConvergence node ? _displayTemporalCriteria(t, metrics, nitrun, format=None, gmres=None, verbose='firstlast', stopAtNan=True) The doc of this function is not yet available, sorry ! Cheers

Luispain commented 2 months ago

Hi Stéphanie,

Thank you for your answer. I tried what you suggest, but GlobalConvergenceHistory node remains empty:

nit = 100; time = 0
time_step = Internal.getNodeFromName(t, 'time_step')
time_step = Internal.getValue(time_step)
for it in range(nit):
    FastS._compute(t, metrics, it, tc)
    if it%modulo_verif == 0:
        FastS._displayTemporalCriteria(t, metrics, it, format='store', gmres=None, verbose='firstlast', stopAtNan=True)
    time += time_step
vincentcasseau commented 1 month ago

Hi Luis,

Once Ivan commits his edits, I'll correct the function that computes the global residuals. For the time being, please define it in your script as

def _calc_global_convergence(t):
    var_list_L2 =['RSD_L2','RSD_L2_diff']
    var_list_Loo=['RSD_oo','RSD_oo_diff']
    neq=5
    if Internal.getValue(Internal.getNodeFromType(t, 'GoverningEquations_t'))=='NSTurbulent':
        neq = 6

    for b in Internal.getBases(t):
        nzones= len(Internal.getZones(b))
        c     = Internal.getNodeByName(b,'GlobalConvergenceHistory')

        zones                 = Internal.getZones(b)
        zone_cong_history     = Internal.getNodeByName(zones[0],'ZoneConvergenceHistory')
        RSD_                  = Internal.getNodeByName(zone_cong_history,'IterationNumber')[1]
        nrec                  = Internal.getValue(zone_cong_history)

        ##NOTE: the value of the ZoneConvergenceHistory node has the number of rec taking at time t
        ##      the shape of the array of IterationNumber is total number of nrec.
        Internal.setValue(c, nrec);

        ##IterationNumber Node
        Internal.createChild(c, 'IterationNumber', 'DataArray_t', numpy.zeros((nrec), Internal.E_NpyInt))
        RSD_b = Internal.getNodeByName(c, 'IterationNumber')[1]
        RSD_b[:] = RSD_[:]

        ##The rest of the nodes in the convergence history
        for var in var_list_L2:
            tmp = numpy.zeros((nrec*neq), numpy.float64)
            Internal.createChild(c, var , 'DataArray_t', tmp) 
        for var in var_list_Loo:
            tmp = numpy.full((nrec*neq), -numpy.inf, dtype=numpy.float64)
            Internal.createChild(c, var , 'DataArray_t', tmp)                    

        total_Ncells = 0
        for z in Internal.getZones(b):
            #if z[0]=='GlobalConvergentHistory': continue
            zone_cong_history     =Internal.getNodeByName(z,'ZoneConvergenceHistory')            
            isCalcCheck           =Internal.getNodeByName(z,'isCalc')

            if isCalcCheck is None:isCalc=1
            else: isCalc=Internal.getValue(isCalcCheck)

            nx = Internal.getValue(z)[0][0]
            ny = Internal.getValue(z)[1][0]
            nz = Internal.getValue(z)[2][0]
            zone_Ncells = nx*ny*nz*isCalc
            total_Ncells += zone_Ncells

            if isCalc==1:
                ##Loo per base
                for var_local in var_list_Loo:
                    RSD_     = Internal.getNodeByName(zone_cong_history,var_local)[1]
                    RSD_b    = Internal.getNodeByName(c,var_local)[1]
                    RSD_b[:] = numpy.maximum(RSD_b, RSD_)

                ##L2 per base (here only sum)
                for var_local in var_list_L2:
                    RSD_     = Internal.getNodeByName(zone_cong_history,var_local)[1]
                    RSD_b    = Internal.getNodeByName(c,var_local)[1]
                    RSD_b[:] += RSD_[:]*zone_Ncells

        ##Average L2 per base
        for var_local in var_list_L2:
            RSD_b    = Internal.getNodeByName(c,var_local)[1]
            RSD_b[:] /= total_Ncells            
    return t

and call it at the end of your script

_calc_global_convergence(t)

which will soon be replaced by

FastS._calc_global_convergence(t)

NB: Stéphanie was right, leave format='store' in your call to FastS._displayTemporalCriteria

Thanks

antjost commented 1 month ago

Hello Luis,

What I might also suggest:

tconv=FastS.create_add_t_converg_hist(t) # creates a PyTree of the convergence histories only. FastS._calc_global_convergence(tconv) # calcs global convergence histories FastS._extractConvergenceHistory(tconv, filename, perZones=True/False, perBases=True/Fase) # exports in Tecplot format

such that one can (1) just load the residuals in tecplot (making their visualization easier) and (2) one can place these functions in the iteration loop to see the residuals on the fly (if it is needed). This is the approach I use. Cheers Antoine

vincentcasseau commented 3 weeks ago

Hi Luis, Could you please confirm that it works now? Commit 62a50fa contains both @IvanMary69's fix to _createConvergenceHistory and a patch for globalConvergenceHistory. Thanks

Luispain commented 3 weeks ago

Hi Vincent,

The new version works great for zonal residuals extraction. However, I still do not get t = FastS.calc_global_convergence(t) working during the simulation (we need the global residuals during simulation).

The case reproducing the issue, adapted from #13 :

# - convergenceHistory (pyTree) -
import numpy
import Converter.PyTree as C
import Generator.PyTree as G
import FastS.PyTree as FastS
import FastC.PyTree as FastC
import Connector.PyTree as X
import Converter.Internal as I
import Initiator.PyTree as Init

npts = 9
dx = 0.5
z = G.cart((0.0,0.0,0.0), (dx,dx,dx), (npts,npts,npts))
C._addBC2Zone(z, 'WALL', 'FamilySpecified:WALL', 'imin')
C._fillEmptyBCWith(z, 'FARFIELD', 'FamilySpecified:FARFIELD', dim=3)
C._addState(z, 'GoverningEquations', 'NSTurbulent')

Init._initConst(z, MInf=0.4, loc='centers')
C._addState(z, MInf=0.4)
t = C.newPyTree(['Base', z])
C._tagWithFamily(t,'FARFIELD')
C._tagWithFamily(t,'WALL')
C._addFamily2Base(t, 'FARFIELD', bndType='BCFarfield')
C._addFamily2Base(t, 'WALL', bndType='BCWall')

import Dist2Walls.PyTree as DTW
walls = C.extractBCOfType(t, 'BCWall')
DTW._distance2Walls(t, walls, loc='centers', type='ortho')

numb = { 'temporal_scheme': 'implicit', 'ss_iteration':3, 'modulo_verif':1}
numz = { 'scheme':'roe', 'slope':'minmod',
    'time_step':0.0007,'time_step_nature':'local', 'cfl':4}
FastC._setNum2Zones(t, numz); FastC._setNum2Base(t, numb)

(t, tc, metrics) = FastS.warmup(t, None)

nit = 10; time = 0

#To remove old ConvergenceHistory nodes 
t = C.rmNodes(t, "ZoneConvergenceHistory")

FastS._createConvergenceHistory(t, nit)

time_step = I.getNodeFromName(t, 'time_step')
time_step = I.getValue(time_step)
for it in range(nit):
    FastS._compute(t, metrics, it, tc)
    FastS.display_temporal_criteria(t, metrics, it, format='store')
    t = FastS.calc_global_convergence(t)
    time += time_step

# time stamp
I.createUniqueChild(t, 'Iteration', 'DataArray_t', value=nit)
I.createUniqueChild(t, 'Time', 'DataArray_t', value=time)
C.convertPyTree2File(t, 'out.cgns')

RSD_L2 = I.getNodeFromName(t,'RSD_L2')[1]
print(f"{RSD_L2=}")

ItNumber = I.getNodeFromName(t,'IterationNumber')[1]
print(f"{ItNumber=}")

Which provokes the error:

Traceback (most recent call last):
  File "fast13.py", line 82, in <module>
    t = FastS.calc_global_convergence(t)
  File "/stck/cassiope/git/Cassiopee/Dist/bin/ld/lib/python3.8/site-packages/FastS/PyTree.py", line 3888, in calc_global_convergence
    RSD_b[:] = RSD_[:]
ValueError: could not broadcast input array from shape (10,) into shape (1,)

Thank you for your support

Best regards Luis