sfepy/sfepy — Issue #38: term tests with manufactured solutions
Author: rc
Created: 2012-01-30 10:24:31+00:00
State: open

Test terms and equations using the method of manufactured solutions (MMS).

Briefly, the MMS involves the following steps:
1. assume a particular solution (e.g. u = sin(x) * cos(y))
2. analytically compute the required right-hand side (source terms) and
boundary conditions for a particular domain
3. compare the numerical solution with the analytical one (possibly for
various mesh refinements).

References (via http://doi.org/):
doi:10.1016/j.jcp.2004.10.036
doi:10.1115/1.1436090

Example:
http://www.imechanica.org/node/1357

How (quite complex!):
- try to use sympy to have an automatic derivation of the source terms.
- annotate somehow the term classes with what they do symbolically using
sympy syntax

First, we should try a simple test case with the Laplace term.

Migrated from http://code.google.com/p/sfepy/issues/detail?id=34

earlier comments

robert.cimrman said, at 2008-05-07T14:27:12.000Z:

The simple test case works, see tests/test_msm_laplace.py, now let us do the complex
stuff. Even this test showed me how the 'per partes' integration in deriving weak
forms can bite, when one forgets the fluxes...

robert.cimrman said, at 2008-05-09T14:18:56.000Z:

The first prototype using SymPy to generate the manufactured right-hand sides is done
(rev. 195:cca307dfbf73, see
http://hg.sympy.org/sfepy/file/cca307dfbf73/tests/test_msm_laplace.py).

Sample output:
... evaluating terms, "<=" is solution, "=>" is the rhs:
... dw_laplace( material, virtual, state )
... symbolic: c * div( grad( u, [x, y, z] ), [x, y, z] )
... using argument map: {'c': 'material', 'u': 'state'}
... <= sin( 3.0 * x ) * cos( 4.0 * y )
... => -25*cos(4*y)*sin(3*x)
... <= ((x - 0.5)**3) * sin( 5.0 * y )
... => -3*(1 - 2*x)*sin(5*y) + 25*(0.5 - x)**3*sin(5*y)
... <= (x**2) + (y**2)
... => 4

Note that the material parameter 'c' of the Laplace term is actually evaluated in the
same way as in a regular computation -> changing it within the input file
tests/test_msm_laplace.py (see material_1) from 1 to 12 leads to:

... evaluating terms, "<=" is solution, "=>" is the rhs:
... dw_laplace( material, virtual, state )
... symbolic: c * div( grad( u, [x, y, z] ), [x, y, z] )
... using argument map: {'c': 'material', 'u': 'state'}
... <= sin( 3.0 * x ) * cos( 4.0 * y )
... => -3E+2*cos(4*y)*sin(3*x)
... <= ((x - 0.5)**3) * sin( 5.0 * y )
... => -36*(1 - 2*x)*sin(5*y) + 3E+2*(0.5 - x)**3*sin(5*y)
... <= (x**2) + (y**2)
... => 48

It works, but there are some simplifying assumptions, namely a single scalar unknown
variable. I also do not like the explicit listing of variables as in
http://hg.sympy.org/sfepy/file/f304324986ff/sfe/terms/termsLaplace.py

symbolic = {'expression': 'c * div( grad( u, [x, y, z] ), [x, y, z] )',
'map' : {'u' : 'state', 'c' : 'material'}}

Nevertheless, the Pandora's box is open!

robert.cimrman said, at 2008-05-16T10:17:48.000Z:

I will remove the [x, y, z] from 'c * div( grad( u, [x, y, z] ), [x, y, z] )', as a
solution is always a function of space (and time), so x, y, z, and t variables will
always be defined.

Other thing: should scalar fields be somehow distinguished from vector fields to
allow the following definition? linear elasticity, small strain tensor, just a
prototype!:
e = 1/2 * (grad( vec( u ) ) + grad( vec( u ) ).T)
instead of
e_ij = 1/2 * (d_j( u_i ) + d_i( u_j ))
Or not? I am not sure yet.

Anyway, I would rather not reinvent SyFi...

ondrej.certik said, at 2008-05-16T11:52:28.000Z:

I don't mind either way.

ondrej.certik said, at 2009-04-01T07:21:58.000Z:

will have to implement it myself
:)
I hoped you will do it
:)
we need something in sympy, so that we can specify all those weakforms in sympy
and let sympy do the thing
so we need a Nabla operator
we need a dot product (I guess just regular "*"?)
etc.