Multi-physics interface

[Jaakko Leppänen]

NOTE: Radial division was added in the output power distribution in version 2.1.10. The description below is for the new format.

One of the major projects in Serpent 2 development is the implementation of a universal multi-physics interface for coupling Serpent to thermal hydraulics and fuel performance codes. This topic is far too extensive to be covered here, and we'll publish more information as the work proceeds. The general idea is that material densities and temperatures can be passed to Serpent without modifications in the main input file, and the output (power distributions, etc.) is written in a separate file that can be easily read and passed back to the coupled code.

The work started a few versions ago, but in update 2.1.5 the routines are getting ready to be used in actual calculations. There are currently three interface types that can be used for defining material density distributions:

1) Continuous distribution by averaging over point-wise values
2) Piece-wise homogeneous distribution in a mesh
3) User-defined functional distribution

The interface is invoked by setting the file name in the main input:

Code: Select all

ifc <file>

One or more interface files can be defined for each problem, for different distributions in coolant and moderator channels, for example. The interface file is a text file, with contents depending on the interface type. The first two lines are common to all types:

Code: Select all

<type> <mat> <out>
<outf> <nz> <zmin> <zmax> <nr>

Parameter <type> is the interface type (1, 2 or 3), <mat> is the name of the material for which the distributions are used and <out> is a flag for printing output (0 = no, 1 = yes). If the flag is not set, the second line is omitted. If set, parameter <outf> determines the file name to which the output is printed, <nz> is the number of axial zones for pin-wise power distributions and <zmin> and <zmax> are the minimum and maximum axial coordinates of the distribution. Radially the pin is divided in <nr> rings with equal volume.

The remaining parameters for type 1 are:

Code: Select all

<dim> <rad> <exp>
<np>
<x1> <y1> <z1> <dens1> <T1>
<x2> <y2> <z2> <dens2> <T2>
...

where <dim> is the dimension of the distribution (1 = z-dependence only, 2 = x- and y-dependence, 3 = x-, y- and z-dependence). The distribution is defined by <np> points with coordinates <x>, <y> and <z> (if dim = 1, only <z> is given, if dim = 2, only <x> and <y> are given). Parameter <dens> is the material density at the point (negative values for mass density, positive values for atomic density) and <T> is the temperature, which is currently not used. Values between the points are calculated by averaging over the listed values. The way the averaging is done depends on parameter <exp>. If <exp> = 1, the average is based on distance, if <exp> = 2, on square distance, and so on. Parameter <rad> determines an exlcusion radius, beyond which the points are not included in the average. If the exclusion radius is too small, some parts of the region may not be properly covered. If the radius is large, more CPU time is required for calculating the average.

Parameters for type 2 are:

Code: Select all

<mtp> <nx> <xmin> <xmax> <ny> <ymin> <ymax> <nz> <zmin> <zmax>
<dens1> <T1>
<dens2> <T2>
...

where <mtp> is the mesh type (1 = Cartesian, the only one currently allowed). The remaining parameters determine the number of mesh cells and the dimensions. The <nx>*<ny>*<nz> cell values follow (in order: x, y, z).

Parameters for type 3 are:

Code: Select all

<np>
<p1>
<p2>
...

The interface contains <np> arbitrary input parameters, that are passed to subroutine userdf.c, along with local coordinates. The subroutine can be modified to include any kind of functional distribution. The interface file has almost no error tolerance, so be sure that all required values are entered and in the correct order.

If the output flag is set, the code calculates power distributions in pin-type objects that are enclosed within the interface material (e.g. fuel pins inside coolant). The output file contains fission powers in cylindrical regions, one region per line. The parameters in each line are:

Code: Select all

<x0> <y0> <z0> <x1> <y1> <z1> <r0> <r1> <P> <dP>

The first six values give the top and bottom center coordinates of the cylindrical region. The next two parameters are the inner and outer radius of the ring, and the last two values are the integral fission power and the associated relative statistical error.

I don't expect the interface to work yet in very complicated geometries (I've tried it mainly in 1D distributions at fuel assembly level). The formats will most likely change in future updates, as more features are added, and the main reason why I'm posting these instructions is that some users have expressed their interest in trying out these capabilities after presentations given on Serpent 2. So feel free to use it (and feedback and new ideas are always welcome), but don't expect too much yet.

The distributions are not shown in the geometry plotter, but there is a makeshift plotting capability implemented in the mesh plotter type 11 (see related post on new mesh plotter capabilities). The plot that is produced is based on collisions, so there is some averaging and statistics involved, but it can be useful for getting an idea on how the code sees the distributions.

I'll post some examples later.

[Jaakko Leppänen]

The input and output formats were changed in update 2.1.10 (added radial binning for the output power distributions).

[Jaakko Leppänen]

A few simple examples of the different interface types can be found at:

http://virtual.vtt.fi/virtual/montecarlo/examples/IFC/

The calculation case is a single 3.7m long BWR fuel assembly, in which the coolant density changes as function of axial position. Input file "ref" is a reference case, in which the distribution is set up by dividing the channel into a number of material cells, each assigned with a different density. Files "bwr1", "bwr2" and "bwr3" use the three types of interfaces (files "dis1.in", "dis2.in" and "dis3.in") to define the same distribution. The density in each case depends only on the axial coordinate, and adding x- and y-dependence complicates the definition to some extent.

While running the calculations I noticed a few problems in the latest version (2.1.10):

1) In debug mode, the calculation terminates in an error message that results from an index check that is done wrong. Without debug mode the calculation runs just fine.
2) Some integral parameters, at least IMP_KEFF and TOT_RR are incorrect. This is a problem I fixed in version 2.1.8, but then managed to re-introduce in the latest update. Detector output and power distribution printed in the interface output file should be OK.
3) The multi-physics interface does not work with universe symmetries.

Also note that the example cases do not work with versions earlier than 2.1.10, since the number of interface input parameters have been changed.

NOTE: Problems 1 and 2 above were fixed in update 2.1.11.

[Jaakko Leppänen]

Temperature was added in the user subroutine in update 2.1.13 and the source file was renamed userifc.c.

[Jaakko Leppänen]

An unstructured mesh based interface using OpenFOAM mesh and file format was added in update 2.1.17. This interface is used similar to the previous interfaces, only the syntax of the file linked with the ifc card is different:

Code: Select all

7 <mat> 0
<rho0> <T0>
<msh_split> <msh_dim> <s0> <sz1> ... <sz_dim>
<points_file>
<faces_file>
<owner_file>
<neighbour_file>
<density_file>
<temperature_file>

where the first line is similar to the other interfaces, only the type is set to 7. The output routines have not yet been implemented, so the output flag must be set to 0. The following line gives the nominal density and temperature (<rho0> and <T0>).

The third line defines the adaptive search mesh. The first parameter (<msh_split>) gives the limiting number of cells, after which the search mesh cell is split into a new mesh. The second parameter (<msh_dim>) gives the mesh dimension, i.e. the maximum number of times a search mesh cell can be split. This value is followed by an equal number of size parameters, giving the size of the search mesh at each level. For example, entry:

Code: Select all

1000 2 10 5

means that the search mesh consists of two levels. The first level is formed by a 10x10x10 mesh. If the number of OpenFOAM cells falling into any of these cells exceeds 1000, the cell is split into a new 5x5x5 mesh. The idea is that the search mesh is automatically refined in a region where the resolution is high, leading to a speed-up in the cell search routine. The methodology is still very much under development.

The remaining entries define the file paths to the OpenFOAM format mesh files: points, faces, owner, neighbour, density and temperature. The last two can be replaced by "-1", in which case the values are replaced by the nominal values. The density file should contain the density factors, i.e. the relative densities compared to the nominal value, and the temperature file the absolute temperatures.

The methodology has been tested in some relatively simple geometries, and a working example (simplified MSR core) can be found at:

http://virtual.vtt.fi/virtual/montecarl ... Interface/

A few things to note:

1) The output routines have not yet been implemented, so the interface can only pass density and temperature data into Serpent.

2) The generation of the search mesh becomes a major bottleneck if the number of dimensions is large (> 2). This is due to a more general problem related to memory allocation in OpenMP parallel mode, and I'll try to do something about it in the next update. Also, the speed-up from the adaptivity is not as significant as expected, so there's a chance that the routine is not working properly.

3) If the OpenFOAM mesh is not of tetrahedral type, it is converted to tet-mesh by splitting the cells into a number of tetrahedrons. So far this conversion procedure has been tested with only a single input file.

4) You may get some warning messages about possibly incorrect file formats. If the calculation runs, these messages can be ingored.

This (and other unstructured mesh based interfaces) are among the most important topics in Serpent 2 development at the moment, so I'm looking forward to any feedback from the users!

[Jaakko Leppänen]

The generation of search mesh should be somewhat faster in version 2.1.18 compared to 2.1.17, due to changes in memory allocation routines.

[vitors]

Dear Jaakko,

Could you please elaborate a little bit more on the "adaptive search mesh"?

I'm sure I'm getting it wrong, but this is what I expected to happen when I pass an OpenFOAM mesh to Serpent:

- Serpent gets the non-structured cells and use these cells to perform its calculations. This is possible due to the information (faces, neighbours, points, etc)
provided about the geometry/mesh.

So, my first question: there is no information on boundaries (the boundary file used in OpenFOAM is missing), what makes me think that for now it is only possible
to provide to Serpent one region mesh and this mesh enclosing surface must match one of the defined surfaces inside Serpent. (Sorry, I was unable to read the mesh information provided in msr1 example in paraview in order to clarify this question).

The second question (probably naive) is: what is the point in using an adaptive mesh in neutronics calculations if the fluid simulation code is unable to work at the same
level of details? (Supposing OpenFOAM is using a fixed mesh)

I just tried to change the msr1 example adding a sourrounding water layer in order to provide two openfoam meshes to Serpent using two ifc directives. I don't even now if it is possible, but I got the following error:

Code: Select all

Fatal error in function PotCorr:

Nuclide hwj3.11t is not flagged for ETTM or DBRC

Simulation aborted.

which doesn't know if is really related to my tests.

Sorry if my questions are to basic (or stupid), but I'm justing starting to learn about Monte Carlo simulations coming from the OpenFOAM side.

My best regards,

Vitor

[vitors]

Hello all,

I got rid of the error mentioned in my first post simply ignoring moderation from water. Despite not physically correct, I am trying to test the ifc interface at this moment.
I'll post again when/if I have useful information.

Sorry for the self-reply,

Vitor

[Jaakko Leppänen]

The interface does not define the cell boundaries, only the density and temperature distributions inside the cells (to be precise, the boundary information is also there -- all faces without neighbour information belong to the boundary, Serpent just doesn't do anything with this information). The actual Monte Carlo geometry, including the cell boundaries, are defined in the Serpent input file, and the temperature and density distributions sit on top of the homogeneous material regions. Essentially the idea is to separate the geometry information from the state-point information, which for the user means that there is no need to subdivide the geometry into separate regions to represent the distributions. The next update will include output for power distributions, tallied using the same mesh, so this information can be passed back into OpenFOAM.

Currently the main limitation in the methodology is that the temperatures of S(a,b) thermal scattering laws cannot be adjusted, which is why you got the error message. So basically the results are physically wrong, but there is so much work to be done that we wanted to get the routines into distribution for testing, while we keep working on the physics model.

The adaptivity works by gradually refining the search mesh in regions where the density of OpenFOAM mesh cells is large. The first parameter in the search mesh input gives the maximum number of OpenFOAM cells inside a search mesh cell, after which the cell is subdivided into a new mesh. So if the input is:

Code: Select all

1000 3 10 5 2

the code first creates a 10x10x10 search mesh, and starts adding the OpenFOAM cells into the structure. Once the number of OpenFOAM cells in a search mesh cell exceeds 1000, the cell is subdivided into a new 5x5x5 search mesh. If one of those cells gets more than 1000 OpenFOAM cells, it is again split into a 2x2x2 mesh. This may not be the optimal approach, because in my opinion it requires too much input from the user, but I think it serves well for testing purposes until we come up with a smarter solution.

I will include more examples after the next update. We are also working on two unstructured mesh-based geometry models, in which also the geometry cell (material region) information is included in the structure. The first model is similarly based on the OpenFOAM mesh, and the second one reads geometry cell information from stereolithography (STL) files that can be written by CAD software.

[orca.blu]

Hi Vitor,

In order to visualize the case with paraview or paraFoam,
you need also the “boundary” file, which is not required by Serpent.

Here is the boundary file for the MSR example:

Code: Select all

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

1
(
    Wall
    {
        type            patch;
        physicalType    patch;
        nFaces          2276;
        startFace       19446;
    }
)

// ************************************************************************* //

You have to arrange the files in the correct folders, e.g.:

costant
>> polyMesh
>> >> boundary
>> >> faces
>> >> neighbour
>> >> owner
>> >> points

0
>> T
>> rhok

You also need the “system” folder and a valid controlDict,
you can probably skip fvSchemes and fvSolution.

Ciao