TreatmentHOSol2Output¶

Description¶

Interpolate a JAGUAR high-order solution onto a user chosen number of output points.

Construction¶

import antares
myt = antares.Treatment('')

Parameters¶

base: Base
The input base.
mesh: Base
The input mesh.

output_npts_dir: int,
output_npts_dir.

Preconditions¶

Postconditions¶

Example¶

import antares
myt = antares.Treatment('')
 myt.execute()

Main functions¶

class antares.treatment.codespecific.jaguar.TreatmentHOSol2Output.TreatmentHOSol2Output¶

execute()¶

Returns:

Return type:

Base

Example¶

"""
This example shows how to process a jaguar solution file.

Read the jaguar file.
Read the Gmsh file.
Create a uniform mesh.
"""
import os
import antares

if not os.path.isdir('OUTPUT'):
    os.makedirs('OUTPUT')

r = antares.Reader('jag')
r['filename'] = os.path.join('..', 'data', 'HIGHORDER', 'covo2D.jag')
base = r.read()

print(base)
print(base[0])
print(base[0][0])

reader = antares.Reader('fmt_gmsh')
reader['filename'] = os.path.join('..', 'data', 'HIGHORDER', 'covo2D.msh')
mesh_real = reader.read()
# mesh_real.attrs['max_vpc'] = 8

t = antares.Treatment('HOSol2Output')
t['base'] = base
t['mesh'] = mesh_real
t['output_npts_dir'] = 4
result = t.execute()

writer = antares.Writer('bin_tp')
writer['base'] = result
writer['filename'] = os.path.join('OUTPUT', 'ex_jaguar_covo2d.plt')
writer.dump()

TreatmentJaguarInit¶

Description¶

Interpolate coordinates of a GMSH grid to a basis based on the jaguar input file p parameter. Then compute a solution with the new coordinates and save it as h5py solution file to be read with the coprocessing version of jaguar.

Construction¶

import antares
myt = antares.Treatment('jaguarinit')

Parameters¶

base: Base
The input base coming from the GMSH reader.

jag_config: antares.utils.high_order.Jaguar
The Jaguar object built from the jaguar input file.
init_func: Callable[[ antares.utils.high_order.HighOrderSol ], list(str)]
A function computing initial variables from coordinates and input parameter. Return a list of variable names. Check the example below.

Example¶

The following example shows how to create an initial solution for the covo/3D example.

import numpy as np
import antares
from antares.utils.high_order import Jaguar

def densityVortexInitialization2D(init_sol):

    with open("2DVortex.dat",'r') as fin:
        lines = [float(i.split('!')[0]) for i in fin.readlines()]

    G   = lines[0]
    Rc  = lines[1]
    xC  = lines[2]
    yC  = lines[3]
    cut = lines[4]

    radiusCut = cut*Rc

    gamma = init_sol.gas['gamma']
    Rgas  = init_sol.gas['Rgas']
    Cp    = init_sol.gas['Cp']

    roInf = init_sol.infState['rho']
    uInf = init_sol.infState['u']
    vInf = init_sol.infState['v']

    roEInf   = init_sol.infState['rhoe']
    presInf  = (roEInf - 0.5 * roInf * (uInf*uInf + vInf*vInf)) * (gamma - 1.0)
    TInf     = presInf / (roInf * Rgas)

    X = init_sol.base['z0'].shared['x']
    Y = init_sol.base['z0'].shared['y']

    R2 = (X-xC)**2 + (Y-yC)**2


    # limit the vortex inside a cut*Rc box.
    expon = np.exp(-0.5*R2/Rc**2)

    u = uInf - (Y-yC) / Rc * G * expon
    v = vInf + (X-xC) / Rc * G * expon

    temp = TInf - 0.5 * (G * G) * np.exp(-R2 / Rc**2) / Cp
    ro   = roInf * (temp / TInf)**( 1. / (gamma - 1.) )
    pres = ro * Rgas * temp
    roE = pres / (gamma - 1.) + 0.5 * ( u*u + v*v ) * ro

    # Build solution array as [ncells, nvars, nSolPts] array.
    vlist = ['rho','rhou','rhov','rhow','rhoe']
    init_sol.base['z0']['0'][vlist[0]]  = np.where(np.sqrt(R2)<=radiusCut, ro, roInf)
    init_sol.base['z0']['0'][vlist[1]] = np.where(np.sqrt(R2)<=radiusCut, ro*u, roInf*uInf)
    init_sol.base['z0']['0'][vlist[2]] = np.where(np.sqrt(R2)<=radiusCut, ro*v, roInf*vInf)
    init_sol.base['z0']['0'][vlist[3]] = 0.
    init_sol.base['z0']['0'][vlist[4]] = np.where(np.sqrt(R2)<=radiusCut, roE, roEInf)

    return vlist


def TGVinit(init_sol):
    gamma = init_sol.gas['gamma']
    Rgas  = init_sol.gas['Rgas']
    Cp    = init_sol.gas['Cp']

    roInf = 1.
    TInf = 351.6
    presInf = 101325.
    u0 = 1
    L = 2*np.pi

    X = init_sol.base['z0'].shared['x']
    Y = init_sol.base['z0'].shared['y']
    Z = init_sol.base['z0'].shared['z']

    u =  u0 * np.sin(X)*np.cos(Y)*np.cos(Z)
    v = -u0 * np.cos(X)*np.sin(Y)*np.cos(Z)
    w = 0.0

    roE = presInf / (gamma - 1.) + 0.5 * ( u*u + v*v + w*w) * roInf

    # Build solution array as [ncells, nvars, nSolPts] array.
    vlist = ['rho','rhou','rhov','rhow','rhoe']
    init_sol.base['z0']['0'][vlist[0]] = roInf
    init_sol.base['z0']['0'][vlist[1]] = roInf*u
    init_sol.base['z0']['0'][vlist[2]] = roInf*v
    init_sol.base['z0']['0'][vlist[3]] = 0.0
    init_sol.base['z0']['0'][vlist[4]] = roE

    return vlist


###################################################"

def main():

    jaguar = Jaguar("input.txt")

    # READER
    # ------

    # Read the GMSH mesh
    r = antares.Reader('fmt_gmsh')
    r['filename'] = jaguar.config['GENERALITIES']['Grid file']
    base_in = r.read()
    base_in.attrs['filename'] = r['filename']


    # TREATMENT
    # ---------

    treatment = antares.Treatment('jaguarinit')

    treatment['base'] = base_in
    treatment['jag_config'] = jaguar
    treatment['init_func'] = TGVinit

    base_out = treatment.execute()


    # WRITER
    # ------

    w = antares.Writer('hdf_jaguar_restart')
    w['base'] = base_out
    w['strategy'] = 'monoproc'
    w.dump()


if __name__ == "__main__":
    main()

Main functions¶

class antares.treatment.codespecific.jaguar.TreatmentJaguarInit.TreatmentJaguarInit¶

execute()¶

Create an initial solution for Jaguar.

Returns:: the unstructured Base obtained from the initial function applied on the input Base coordinates.
Return type:: Base

Example¶

"""
This example creates a (restart) initial solution for a jaguar computation with coprocessing.
It is based on the covo3D jaguar test case.

Execute it directly in the PATH_TO_JAGUAR/tests/covo/3D directory.
Tested on kraken with:
1/ source /softs/venv_python3/vant_ompi401/bin/activate
2/ source <path_to_antares>/antares.env
3/ python create_jaguar_init_sol.py

Output should be:
```
===================================================================
                          ANTARES 1.14.0
                    https://cerfacs.fr/antares
Copyright 2012-2019, CERFACS. This software is licensed.
===================================================================

#### CREATE JAGUAR INSTANCE... ####
>>> Read input file input.txt...
>>> Mesh file cube_8_p4.h5 Found.
>>> Light load mesh file cube_8_p4.h5 ...
>>> Mesh file cube_8_p4_vtklag.h5 Found.
>>> Light load mesh file cube_8_p4_vtklag.h5 ...
input x_mesh.shape:: (512, 8)
output x_new.shape:: (512, 125)
> 5 var in sol: ['rho', 'rhou', 'rhov', 'rhow', 'rhoe']
> 0 var in additionals: []
>>> XMF file restart_0000000.xmf written.
>>> Solution file restart_0000000.jag written.
> 5 var in sol: ['rho', 'rhou', 'rhov', 'rhow', 'rhoe']
> 0 var in additionals: []
>>> XMF file restart_0000000_vtklag.xmf written.
>>> Solution file restart_0000000_vtklag.jag written.
```

4/ Read the file restart_0000000.xmf (solution pts) or restart_0000000_vtklag.xmf (output pts) in ParaView.
"""
import os

import numpy as np

import antares
from antares.utils.high_order import Jaguar

def pulse2D(mesh, init_sol):
    ro = init_sol.infState['rho']
    uInf = init_sol.infState['u']
    vInf = init_sol.infState['v']

    p0=1e5
    p0prime=1e3
    x0 = 0.05
    y0 = 0.05

    pres=p0+p0prime*np.exp((-np.log10(2.0)/.0001)*((X-x0)*(X-x0)+(Y-y0)*(Y-y0)))
    roE=pres/(init_sol.infState['gamma']-1)+0.5*ro*(uInf*uInf+vInf*vInf)

    # Build solution array as [ncells, nvars, nSolPts] array.
    vlist = ['rho','rhou','rhov','rhow','rhoe']
    init_sol.base['z0']['0'][vlist[0]] = ro
    init_sol.base['z0']['0'][vlist[1]] = init_sol.infState['rhou']
    init_sol.base['z0']['0'][vlist[2]] = init_sol.infState['rhov']
    init_sol.base['z0']['0'][vlist[3]] = 0.
    init_sol.base['z0']['0'][vlist[4]] = roE

    return vlist

def densityVortexInitialization2D(init_sol):

    with open(os.path.join('..', 'data', 'GMSH', '2DVortex.dat'),'r') as fin:
        lines = [float(i.split('!')[0]) for i in fin.readlines()]

    G   = lines[0]
    Rc  = lines[1]
    xC  = lines[2]
    yC  = lines[3]
    cut = lines[4]

    radiusCut = cut*Rc

    gamma = init_sol.gas['gamma']
    Rgas  = init_sol.gas['Rgas']
    Cp    = init_sol.gas['Cp']

    roInf = init_sol.infState['rho']
    uInf = init_sol.infState['u']
    vInf = init_sol.infState['v']

    roEInf   = init_sol.infState['rhoe']
    presInf  = (roEInf - 0.5 * roInf * (uInf*uInf + vInf*vInf)) * (gamma - 1.0)
    TInf     = presInf / (roInf * Rgas)

    X = init_sol.base['z0'].shared['x']
    Y = init_sol.base['z0'].shared['y']

    R2 = (X-xC)**2 + (Y-yC)**2


    # limit the vortex inside a cut*Rc box.
    expon = np.exp(-0.5*R2/Rc**2)

    u = uInf - (Y-yC) / Rc * G * expon
    v = vInf + (X-xC) / Rc * G * expon

    temp = TInf - 0.5 * (G * G) * np.exp(-R2 / Rc**2) / Cp
    ro   = roInf * (temp / TInf)**( 1. / (gamma - 1.) )
    pres = ro * Rgas * temp
    roE = pres / (gamma - 1.) + 0.5 * ( u*u + v*v ) * ro

    # Build solution array as [ncells, nvars, nSolPts] array.
    vlist = ['rho','rhou','rhov','rhow','rhoe']
    init_sol.base['z0']['0'][vlist[0]]  = np.where(np.sqrt(R2)<=radiusCut, ro, roInf)
    init_sol.base['z0']['0'][vlist[1]] = np.where(np.sqrt(R2)<=radiusCut, ro*u, roInf*uInf)
    init_sol.base['z0']['0'][vlist[2]] = np.where(np.sqrt(R2)<=radiusCut, ro*v, roInf*vInf)
    init_sol.base['z0']['0'][vlist[3]] = 0.
    init_sol.base['z0']['0'][vlist[4]] = np.where(np.sqrt(R2)<=radiusCut, roE, roEInf)

    return vlist


def TGVinit(init_sol):
    gamma = init_sol.gas['gamma']
    Rgas  = init_sol.gas['Rgas']
    Cp    = init_sol.gas['Cp']

    roInf = 1.
    TInf = 351.6
    presInf = 101325.
    u0 = 1
    L = 2*np.pi

    X = init_sol.base['z0'].shared['x']
    Y = init_sol.base['z0'].shared['y']
    Z = init_sol.base['z0'].shared['z']

    u =  u0 * np.sin(X)*np.cos(Y)*np.cos(Z)
    v = -u0 * np.cos(X)*np.sin(Y)*np.cos(Z)
    w = 0.0

    roE = presInf / (gamma - 1.) + 0.5 * ( u*u + v*v + w*w) * roInf

    # Build solution array as [ncells, nvars, nSolPts] array.
    vlist = ['rho','rhou','rhov','rhow','rhoe']
    init_sol.base['z0']['0'][vlist[0]] = roInf
    init_sol.base['z0']['0'][vlist[1]] = roInf*u
    init_sol.base['z0']['0'][vlist[2]] = roInf*v
    init_sol.base['z0']['0'][vlist[3]] = 0.0
    init_sol.base['z0']['0'][vlist[4]] = roE

    return vlist


###################################################"
def main():

    jaguar_input_file = os.path.join('..', 'data', 'GMSH', 'input.txt')
    jaguar = Jaguar(jaguar_input_file)

    # READER
    # ------
    # Read the GMSH mesh
    r = antares.Reader('fmt_gmsh')
    r['filename'] = os.path.join('..', 'data', 'GMSH', 'cube_8.msh')
    base_in = r.read()
    base_in.attrs['filename'] = r['filename']

    # TREATMENT
    # ---------
    treatment = antares.Treatment('jaguarinit')

    treatment['base'] = base_in
    treatment['jag_config'] = jaguar
    treatment['init_func'] = densityVortexInitialization2D
    # treatment['init_func'] = TGVinit # TGV function may also be tested

    base_out = treatment.execute()

    # WRITER
    # ------
    w = antares.Writer('hdf_jaguar_restart')
    w['base'] = base_out
    w['strategy'] = 'monoproc'
    w.dump()


if __name__ == "__main__":
    main()

High Order Tools¶

class antares.utils.high_order.HighOrderTools¶

LagrangeGaussLegendre(p)¶

cell_sol_interpolation(sol_in, p_in, x_out)¶

Interpolate sol_in, p_in on a mesh based on 1D isoparam location x_out

Arguments:: sol_in {array} – Input solution p_in {array} – array of order of each cell of sol_in x_out {array} – 1D location of new interpolation points in [0,1]

getDerLi1D(x)¶: Compute derivative of Lagrange basis. :x: input 1D points where data is known return: matrix : dL1 .. dLn

x1 [dL1(x1) dLn(x1)] .. [ …. ] xn [dL1(xn) .. dLn(xn)]

getDerLi2D(x)¶

Compute the partial derivative values of the 2D Lagrange basis at its construction points x.

Detail: To match the jaguar decreasing x order we have to flip the Li(x) before the tensor product, ie: dLidx = flip( dLi(x) ) * Li(y) dLidy = flip( Li(x) ) * dLi(y) = Li(x) * dLi(y) because Li = np.eye = identity (np.eye is used instead of getLi1D because this is the result of the Lagrange basis on its contruction points.)

getDerLi3D(x)¶

Compute the partial derivative values of the 3D Lagrange basis at its construction points x.

Detail: To match the jaguar decreasing x order we have to flip the Li(x) before the tensor product, ie: dLidx = flip( dLi(x) ) * Li(y) * Li(z) dLidy = flip( Li(x) ) * dLi(y) * Li(z) = Li(x) * dLi(y) * Li(z) because Li = np.eye = identity dLidz = flip( Li(x) ) * Li(y) * dLi(z) (np.eye is used instead of getLi1D because this is the result of the Lagrange basis on its contruction points.)

getLi1D(x_in, x_out)¶

getLi2D(x_in, x_out, flatten='no')¶: Return a matrix to interpolate from 2D based on 1D locations x_in to 2D based on 1D locations x_out

getLi3D(x_in, x_out, flatten='no')¶: Return a matrix to interpolate from 3D based on 1D locations x_in to 3D based on 1D locations x_out

ndim = 0¶

class antares.utils.high_order.HighOrderMesh(base=None)¶

Class to handle high order mesh with Jaguar.

What you can do: - read a gmsh hexa mesh (order 1 or 2) and write a new linear mesh with a given basis to xdmf format.

add_output_mesh(basis, p_out)¶

create_from_GMSH()¶

Read the GMSH mesh with antares and generate x,y(,z)_mesh arrays of input mesh vertices. Warning: this is still a monoproc function.

Set self.x/y/z_mesh arrays

get_names()¶

class antares.utils.high_order.HighOrderSol(parent)¶

Solution of a high order computation with several variables.

create_from_coprocessing(mesh)¶

Link a list of mesh to the solution and set the jagSP as default

Args:: mesh (HighOrderMesh): HighOrder Mesh instance.

create_from_init_function(mesh, init_func)¶

Apply init_func on mesh coordinates to generate a initial solution field.

Args:: mesh (OutputMesh): Instance of a OutputMesh init_func (function): Function taking in argument a HighOrderSol and returning an ordered list of variables.

set_active_mesh(name)¶

class antares.utils.high_order.OutputMesh(parent)¶

Class to handle high order mesh with Jaguar What you can do: - read a gmsh hexa mesh (order 1 or 2) and write a new linear mesh with a given basis to xdmf format. -

invert(list1, list2)¶: Given two maps (lists) from [0..N] to nodes, finds a permutations between them. :arg list1: a list of nodes. :arg list2: a second list of nodes. :returns: a list of integers, l, such that list1[x] = list2[l[x]] source:https://github.com/firedrakeproject/firedrake/blob/661fc4597000ccb5658deab0fa99c3f7cab65ac3/firedrake/paraview_reordering.py

light_load()¶: Load only some scalar parameter from existing mesh file

set_output_basis(basis, p=None, custom_1D_basis=None)¶: Set the output basis. :basis: - ‘jagSP’ : Use SP as OP - ‘jagSPextra’ : use SP as OP and add extrapolated points at 0 and 1 (in iso coor) to have a continuous render - ‘vtklag’ : same as above but boundaries are extended to [0,1] - ‘custom’ : use custom_1D_basis as interpolation basis. :p: basis order, if None use the attached solution order

vtk_hex_local_to_cart(orders)¶: Produces a list of nodes for VTK’s lagrange hex basis. :arg order: the three orders of the hex basis. :return a list of arrays of floats. source:https://github.com/firedrakeproject/firedrake/blob/661fc4597000ccb5658deab0fa99c3f7cab65ac3/firedrake/paraview_reordering.py

write()¶

write_h5()¶: Write high order mesh to different available format

write_xmf()¶: Write xmf file to read the mesh without solution

class antares.utils.high_order.SolWriter(base=None)¶

write_monoproc()¶

Write a single-file solution in single-proc mode.

Mainly used up to now to write init solution.

write_parallel(restart_mode='h5py_ll')¶

Parallel HDF5 write: multi proc to single file

Args:: restart_mode (str, optional): See below. Defaults to ‘h5py_ll’. basis (str, optional): Output basis. Coprocessing returns solution at solution points, but for visualization purpose ‘vtklag’ may be use to generate uniform basis and so be able to use vtk Lagrange cells. Defaults to ‘jagSP’.

Several restart modes have been tested: - h5py_NOTsorted_highAPI: ‘h5py_hl_ns’ - h5py_sorted_highAPI: ‘h5py_hl’ - h5py_sorted_lowAPI: ‘h5py_ll’ - fortran_sorted_opt2: ‘fortran’ In sorted results, the efficiency is from best to worst: fortran > h5py_ll > h5py hl

Here, only h5py low level API (restart_mode=h5py_ll) is used because it has an efficiency near of the pure fortran, but it is much easier to modify and debug. For final best performance, a pure fortran code can be derived quite easily based of existing one from the python h5py_ll code.

class antares.utils.high_order.Jaguar(input_fname)¶

class CoprocOutput¶

Use the formalism <name> and <name>_size when you add array to be able to use jaguar_to_numpy function

conn¶: Structure/Union member

conn_size¶: Structure/Union member

coor¶: Structure/Union member

coor_size¶: Structure/Union member

data¶: Structure/Union member

data_size¶: Structure/Union member

extra_vars_size¶: Structure/Union member

iter¶: Structure/Union member

rms_vars_size¶: Structure/Union member

time¶: Structure/Union member

vars_names¶: Structure/Union member

vars_size¶: Structure/Union member

class CoprocParam¶

Parameters are relative to each partition

antares_io¶: Structure/Union member

cp_ite¶: Structure/Union member

cp_type¶: Structure/Union member

input_file¶: Structure/Union member

input_file_size¶: Structure/Union member

nTot_OP¶: Structure/Union member

n_OP¶: Structure/Union member

n_OP_dir¶: Structure/Union member

n_extra_vars¶: Structure/Union member

n_rms_vars¶: Structure/Union member

ndim¶: Structure/Union member

need_mesh¶: Structure/Union member

nloc_SP¶: Structure/Union member

nloc_cell¶: Structure/Union member

ntot_SP¶: Structure/Union member

ntot_cell¶: Structure/Union member

nvars¶: Structure/Union member

class CoprocRestart¶

Use the formalism <name> and <name>_size when you add array to be able to use jaguar_to_numpy function

c_l2g¶: Structure/Union member

c_l2g_size¶: Structure/Union member

iter¶: Structure/Union member

maxPolynomialOrder¶: Structure/Union member

nloc_cell¶: Structure/Union member

nprocs¶: Structure/Union member

p¶: Structure/Union member

p_size¶: Structure/Union member

sol¶: Structure/Union member

sol_size¶: Structure/Union member

time¶: Structure/Union member

PostProcessing()¶

PreProcessing()¶

Run()¶

coprocessing(cp_type)¶

debug(double_arr, int_arr, int_scal)¶

Callback function to be used for quick debuging

Args:

arr (ct.POINTER): ctypes pointer passing from fortran, may be in or float size (ct.POINTER(ct.c_int)): ctypes pointer to c_int passing from fortran giving the size of arr

Usage:

1/ Add ‘use coprocessing’ in the fortran file you want to debug 2/ Add a ‘call callback_debug(double_array, int_array, int_scal)’ where you want.

Example of call on fortran side in subroutine computeAverageVars:

integer, dimension(2) :: tmp tmp(1) = Mesh%nTotSolutionPoints tmp(2) = userData%n_RMS_vars call callback_pydebug(RMS%averageVariable, tmp, 2)

Giving the following output for the COVO3D:: Iter Time deltaT 1 0.23119E-05 0.23119E-05 pydebug [64000 2] (64000, 2)

Then do whatever you want on python side (Antares treatment and output, etc.)

Warning: there is no optional argument here, you must pass 3 arguments on the fortran side. So you can:

A/ either on fortran side create a temporary variable to fit the debug function call B/ or on python side modify:

1/ the call of the debug function above 2/ the definition of ptrf2 = ct.CFUNCTYPE(..) in the self.Run() method.

getLi1D(x_in, x_out)¶

getLi2D(x_in, x_out)¶

getLi3D(x_in, x_out, flatten='no')¶: Return a matrix to interpolate from 3D based on 1D locations x_in to 3D based on 1D locations x_out

init_mesh_and_sol(additional_mesh={})¶

Create and initialize HighOrderMesh and HighOrderSol instance for coprocessing.

Args:

additional_mesh (list, optional): List of dictionnary to add optional output mesh basis.: Keys are ‘type’ and ‘p’, if ‘p’ is not created string input file p is used. Defaults to None.

python_str_to_fortran(s)¶: Source: http://degenerateconic.com/fortran-json-python/

txt2ini(fname)¶: Convert txt input file in ini file format to be read by the parser later.