Tutorial

The goal of SOLVCON is to help code developers to focus on the numerical algorithms. These computing code can be written in C or other high-speed languages and interfaced with SOLVCON. SOLVCON has a general work flow that support mesh loaders (Gmsh, FLUENT Gambit (R), and CUBIT), MPI, and VTK. These supportive functionalities are ready to help developing problem solvers.

Set up the Environment

Assume:

  • SOLVCON is compiled without problems. See install for more information.
  • The compiled SOLVCON is located at $SCSRC.
  • You are using bash.

Usually we don’t install SOLVCON into the OS, but use environment variable to enable it, so that it’s easier to modify the source code:

export PYTHONPATH=$SCSRC:$PYTHONPATH

And then the following command:

python -c "import solvcon; print solvcon"

should show you the correct location of SOLVCON.

There are various examples located at $SCSRC/examples. To follow the examples, you need to:

  • Install ParaView. On a Debian/Ubuntu, you can do it by executing:

    sudo apt-get install paraview
  • Obtain example data. You can do it by executing:

    scons --get-scdata
    

    in $SCSRC.

More information of the verification examples can be found in Verification.

Configuration

SOLVCON will find each of the solvcon.ini files from current working directory toward the root, and use their settings. Three settings are recognized in [SOLVCON] section:

  • APPS: Equivelent to the environment variable SOLVCON_APPS.
  • LOGFILE: Equivelent to the environment variable SOLVCON_LOGFILE.
  • PROJECT_DIR: Equivelent to the environment variable SOLVCON_PROJECT_DIR. Can be set to empty, which indicates the path where the configuration file locates.

The configurable environment variables:

  • SOLVCON_PROJECT_DIR: the directory holds the applications.
  • SOLVCON_LOGFILE: filename for solvcon logfile.
  • SOLVCON_APPS: names of the available applications, seperated with semi-colon. There should be no spaces.
  • SOLVCON_FPDTYPE: a string for the numpy dtype object for floating-point. The default fpdtype to be used is float64 (double).
  • SOLVCON_INTDTYPE: a string for the numpy dtype object for integer. The default intdtype to be used is int32.
  • SOLVCON_MPI: flag to use MPI.

Mesh Generation (to Be Added)

Before solving any PDE, you need to define the discretized spatial domain of the problem by generating the mesh.

Problem Solver

To demonstrate how to develop a problem solver, SOLVCON provides a “fake” one in solvcon.parcel.fake. The package implements a trivial and meaningless algorithm which is easy to validate. The related files are all in the directory $SCSRC/solvcon/parcel/fake. You can follow the source code and test cases to learn about how to write a problem solver with SOLVCON.

There are two modules, solver and _algorithm, inside that package. They provides two classes: FakeSolver and FakeAlgorithm. The former is the higher-level API and purely written in Python. The latter is implemented with Cython to call low-level C code. The real number-crunching code is written in C:

void fake_calc_soln(sc_mesh_t *msd, sc_fake_algorithm_t *alg)

(Jump to the code listing). Let u_j^n be the solution at j-th cell and n-th time step, and v_j be the volume at j-th cell. This function advances each value in the solution array (sc_fake_algorithm_t.sol and sc_fake_algorithm_t.soln) by using the following expression:

u_j^n = u_j^{n-\frac{1}{2}} + \frac{\Delta t}{2}v_j

Total number of values per cell is set in sc_fake_algorithm_t.neq.

void fake_calc_dsoln(sc_mesh_t *msd, sc_fake_algorithm_t *alg)

(Jump to the code listing). Let (u_{x_{\mu}})_j^n be the x_{\mu} component of the gradient of u_j^n, and (c_{\mu})_j be the x_{\mu} component of the centroid of the j-th cell. \mu = 1, 2 or \mu = 1,
2, 3. This function advances each value in the solution gradient array (sc_fake_algorithm_t.dsol and sc_fake_algorithm_t.dsoln) by using the following expression:

(u_{x_{\mu}})_j^n =
  (u_{x_{\mu}})j^{n-\frac{1}{2}} + \frac{\Delta t}{2}(c_{\mu})_j

Total number of values per cell is set in sc_fake_algorithm_t.neq.

The Python/Cython/C hybrid style may seem complicated, but it is important for performance. The two C functions are wrapped with the Cython methods FakeAlgorithm.calc_soln and FakeAlgorithm.calc_dsoln, respectively. Then, the higher level FakeSolver will use the lower-level FakeAlgorithm to connect the underneath numerical algorithm to the supportive functionalities prepared in SOLVCON.

fake.solver

This is the higher level module implemented in Python.

class solvcon.parcel.fake.solver.FakeSolver(blk, neq=None, **kw)

This class represents the Python side of a demonstration-only numerical method. It instantiates a FakeAlgorithm object. Computation-intensive tasks are delegated to the algorithm object.

Inheritance diagram of FakeSolver

__init__(blk, neq=None, **kw)

Constructor of FakeSolver.

A Block is a prerequisite:

>>> from solvcon.testing import create_trivial_2d_blk
>>> blk = create_trivial_2d_blk()

But the constructor also needs to know how many equations (i.e., number of variables per cell). We must provide the neq argument:

>>> _ = FakeSolver(blk)
Traceback (most recent call last):
    ...
ValueError: neq must be int (but got None)
>>> _ = FakeSolver(blk, neq=1.2)
Traceback (most recent call last):
    ...
ValueError: neq must be int (but got 1.2)
>>> svr = FakeSolver(blk, neq=1)

Each time step is composed by two sub time steps, as the CESE method requires:

>>> svr.substep_run
2

The constructor will create four solution arrays (without initialization):

>>> [(type(getattr(svr, key)), getattr(svr, key).shape)
...     for key in ('sol', 'soln', 'dsol', 'dsoln')]
...     
[(<type 'numpy.ndarray'>, (6, 1)), (<type 'numpy.ndarray'>, (6, 1)),
(<type 'numpy.ndarray'>, (6, 1, 2)), (<type 'numpy.ndarray'>, (6, 1,
2))]
neq = None

Number of equations, or number of variables on each cell. Must be an int.

create_alg()

Create a FakeAlgorithm object.

>>> from solvcon.testing import create_trivial_2d_blk
>>> blk = create_trivial_2d_blk()
>>> svr = FakeSolver(blk, neq=1)
>>> isinstance(svr.create_alg(), _algorithm.FakeAlgorithm)
True
_MMNAMES = ['update', 'calcsoln', 'ibcsoln', 'calccfl', 'calcdsoln', 'ibcdsoln']

See solvcon.solver.MeshSolver._MMNAMES.

The following six methods build up the numerical algorithm. They are registered into _MMNAMES with the present order.

update(worker=None)

Update the present solution arrays (sol and dsol) with the contents in the next solution arrays (dsol and dsoln).

>>> from solvcon.testing import create_trivial_2d_blk
>>> blk = create_trivial_2d_blk()
>>> svr = FakeSolver(blk, neq=1)
>>> # initialize with different solution arrays.
>>> svr.sol.fill(0)
>>> svr.soln.fill(2)
>>> svr.dsol.fill(0)
>>> svr.dsoln.fill(2)
>>> (svr.sol != svr.soln).all()
True
>>> (svr.dsol != svr.dsoln).all()
True
>>> # update and then solution arrays become the same.
>>> svr.update()
>>> (svr.sol == svr.soln).all()
True
>>> (svr.dsol == svr.dsoln).all()
True
calcsoln(worker=None)

Advance sol to soln. The calculation is delegated to FakeAlgorithm.calc_soln.

>>> # build a block before creating a solver.
>>> from solvcon.testing import create_trivial_2d_blk
>>> blk = create_trivial_2d_blk()
>>> # create a solver.
>>> svr = FakeSolver(blk, neq=1)
>>> # initialize the solver.
>>> svr.sol.fill(0)
>>> svr.soln.fill(0)
>>> svr.dsol.fill(0)
>>> svr.dsoln.fill(0)
>>> # run the solver.
>>> ret = svr.march(0.0, 0.01, 100)
>>> # calculate and compare the results in soln.
>>> soln = svr.soln[svr.blk.ngstcell:,:]
>>> clvol = np.empty_like(soln)
>>> clvol.fill(0)
>>> for iistep in range(200):
...     clvol[:,0] += svr.blk.clvol*svr.time_increment/2
>>> # compare.
>>> (soln==clvol).all()
True
ibcsoln(worker=None)

Interchange BC for the soln array. Only used for parallel computing.

calccfl(worker=None)

Calculate the CFL number. For FakeSolver, this method does nothing.

calcdsoln(worker=None)

Advance dsol to dsoln. The calculation is delegated to FakeAlgorithm.calc_dsoln.

>>> # build a block before creating a solver.
>>> from solvcon.testing import create_trivial_2d_blk
>>> blk = create_trivial_2d_blk()
>>> # create a solver.
>>> svr = FakeSolver(blk, neq=1)
>>> # initialize the solver.
>>> svr.sol.fill(0)
>>> svr.soln.fill(0)
>>> svr.dsol.fill(0)
>>> svr.dsoln.fill(0)
>>> # run the solver.
>>> ret = svr.march(0.0, 0.01, 100)
>>> # calculate and compare the results in dsoln.
>>> dsoln = svr.dsoln[svr.blk.ngstcell:,0,:]
>>> clcnd = np.empty_like(dsoln)
>>> clcnd.fill(0)
>>> for iistep in range(200):
...     clcnd += svr.blk.clcnd*svr.time_increment/2
>>> # compare.
>>> (dsoln==clcnd).all()
True
ibcdsoln(worker=None)

Interchange BC for the dsoln array. Only used for parallel computing.

FakeSolver defines the following solution arrays:

FakeSolver.sol = None

This is the “present” solution array (numpy.ndarray) for the algorithm with two sub time step.

FakeSolver.soln = None

This is the “next” or “new” solution array (numpy.ndarray) for the algorithm with two sub time step.

FakeSolver.dsol = None

This is the “present” solution gradient array (numpy.ndarray) for the algorithm with two sub time step.

FakeSolver.dsoln = None

This is the “next” or “new” solution gradient array (numpy.ndarray) for the algorithm with two sub time step.

fake._algorithm

This is the lower level module implemented in Cython. It is composed by two files. $SCSRC/solvcon/parcel/fake/_algorithm.pxd declares the data structure for C. $SCSRC/solvcon/parcel/fake/_algorithm.pyx defines the wrapping code.

C API Declaration

The Cython file _algorithm.pxd defines the following type for the low-level C functions to access the data in a FakeAlgorithm that proxies FakeSolver:

sc_fake_algorithm_t
int neq

This should be set to FakeSolver.neq.

double time

This should be set to MeshSolver.time.

double time_increment

This should be set to MeshSolver.time_increment.

double* sol

This should point to the 0-th cell of FakeSolver.sol. Therefore the address of ghost cells is smaller than sc_fake_algorithm_t.sol.

double* soln

This should point to the 0-th cell of FakeSolver.soln. Therefore the address of ghost cells is smaller than sc_fake_algorithm_t.soln.

double* dsol

This should point to the 0-th cell of FakeSolver.dsol. Therefore the address of ghost cells is smaller than sc_fake_algorithm_t.dsol.

double* dsoln

This should point to the 0-th cell of FakeSolver.dsoln. Therefore the address of ghost cells is smaller than sc_fake_algorithm_t.dsoln.

Wrapper Class

class solvcon.parcel.fake._algorithm.FakeAlgorithm

This class wraps around the C portion of the numerical method.

setup_algorithm(svr)

A FakeAlgorithm object shouldn’t allocate memory. Instead, a FakeSolver object should allocate the memory and pass the solver into the algorithm.

calc_soln()

Wraps the C functions fake_calc_soln() (where the algorithm is defined).

calc_dsoln()

Wraps the C functions fake_calc_dsoln() (where the algorithm is defined).

C Code Listing

solvcon/parcel/fake/src/fake_calc_soln.c

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
/*
 * Copyright (c) 2008, Yung-Yu Chen <yyc@solvcon.net>
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * - Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimer.
 * - Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * - Neither the name of the SOLVCON nor the names of its contributors may be
 *   used to endorse or promote products derived from this software without
 *   specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include <Python.h>

#include "mesh.h"
#include "_algorithm.h"

void sc_fake_calc_soln(sc_mesh_t *msd, sc_fake_algorithm_t *alg) {
    double *psol, *psoln, *pclvol;
    int icl, ieq;
    psol = alg->sol;
    psoln = alg->soln;
    pclvol = msd->clvol;
    for (icl=0; icl<msd->ncell; icl++) {
        for (ieq=0; ieq<alg->neq; ieq++) {
            psoln[ieq] = psol[ieq] + pclvol[0] * alg->time_increment / 2.0;
        };
        psol += alg->neq;
        psoln += alg->neq;
        pclvol += 1;
    };
};

// vim: set ts=4 et:

solvcon/parcel/fake/src/fake_calc_dsoln.c

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
/*
 * Copyright (c) 2008, Yung-Yu Chen <yyc@solvcon.net>
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * - Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimer.
 * - Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * - Neither the name of the SOLVCON nor the names of its contributors may be
 *   used to endorse or promote products derived from this software without
 *   specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include <Python.h>

#include "mesh.h"
#include "_algorithm.h"

void sc_fake_calc_dsoln(sc_mesh_t *msd, sc_fake_algorithm_t *alg) {
    double *pdsol, *pdsoln, *pclcnd;
    int icl, ieq, idm;
    pdsol = alg->dsol;
    pdsoln = alg->dsoln;
    pclcnd = msd->clcnd;
    for (icl=0; icl<msd->ncell; icl++) {
        for (ieq=0; ieq<alg->neq; ieq++) {
            for (idm=0; idm<msd->ndim; idm++) {
                pdsoln[idm] = pdsol[idm]
                            + pclcnd[idm] * alg->time_increment / 2.0;
            };
            pdsol += msd->ndim;
            pdsoln += msd->ndim;
        };
        pclcnd += msd->ndim;
    };
};

// vim: set ts=4 et: