gridlod
is a python code for computing localized element correctors
used in the localized orthogonal decomposition (LOD) method for
numerical homogenization. The code works for any number of dimensions
d (also greater than 3) but is limited to computations on the
hypercube [0,1]^d on a structured grid as coarse mesh, and another
structured grid (refinement of the former) as fine mesh.
The code is written with consideration of the memory usage. Thus, the Petrov--Galerkin formulation is used. Also, fine-scale quantities and coarse-scale quantities are distinguished to be able to discard as much fine-scale information as early as possible to reduce the memory usage.
These packages might be needed:
scikit-sparse
numpy
scipy
You can find some examples on how to use it in the test directory
test
. See e.g. the example in gridlod/test/test_pgexamples.py
. The
example is taken from Daniel Peterseim's paper, Variational
Multiscale Stabilization and the Exponential Decay of Correctors.
The code has been developed in stages, and the terminology is not consistent. A few pointers follow.
world
means the unit hypercube in question.patch
means a subdomain that is also a hypercube but where one particular element is special (typically the center element of the patch for corrector problems).
Hypercubes are described by arrays. E.g. NPatch = np.array([5,2,6])
is a 5 times 2 time 6 element hypercube.
To allow for two meshes, we talk about Coarse and Fine
quantities. NWorldCoarse
for example is the extent of the world in
coarse elements, and NWorldFine
the extent in fine elements. The
elementwise ratio NWorldFine/NWorldCoarse
is consistently called
NCoarseElement
and should be all integers.
All domains are hypercubes and can thus generally be described by its
extent (typically NPatchCoarse
) and a starting index, cartesian
(iPatchCoarse
) or lexicographical (with varying names).
When data (e.g. elementwise coefficient or nodewise solution vector)
is stored linearly, it is stored with the first index running
fastest. E.g.: aFine
is typically the coefficient. aFine[0]
is
coefficient at fine index (0, 0, 0), aFine[1]
at (1, 0 0), and
so on.
There is a lot of index arithmetic going on to be able to do bulk
operations with the numpy routines. The util.py
file contains
several routines for constructing subpatches or finding a
lexicographical index from a cartesian coordinate and so on. The basis
of many routines is pIndexMap
.
The base consists of these files:
fem.py
contains code for the assembly of finite element matrices.interp.py
contains code for interpolation operators.world.py
contains the World and Patch class.coef.py
contains code for localizing a coefficient.femsolver.py
contains code for solving a reference FEM problem.linalg.py
contains code for solving linear systems.util.py
contains code for manipulating grid indices.lod.py
contains code for solving corrector problems.pglod.py
contains code for assembling PG-LOD matricesfunc.py
contains code for describing Q0 and Q1 functions.transport.py
contains code for computing fluxes (not well-maintained)
The test and integration directories contain some tests that can be
illustrative. Use e.g. nosetests
to run the tests.
- Fredrik Hellman
- Tim Keil