term tests with manufactured solutions

[GoogleCodeExporter]

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.

Original issue reported on code.google.com by robert.c...@gmail.com on 28 Apr 2008 at 2:34

[GoogleCodeExporter]

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...

Original comment by robert.c...@gmail.com on 7 May 2008 at 2:27

• Changed state: Started

[GoogleCodeExporter]

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!

Original comment by robert.c...@gmail.com on 9 May 2008 at 2:18

[GoogleCodeExporter]

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...

Original comment by robert.c...@gmail.com on 16 May 2008 at 10:17

[GoogleCodeExporter]

I don't mind either way.

Original comment by ondrej.c...@gmail.com on 16 May 2008 at 11:52

[GoogleCodeExporter]

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

Original comment by ondrej.c...@gmail.com on 1 Apr 2009 at 7:21

[GoogleCodeExporter]

Migrated to http://github.com/sfepy/sfepy/issues/38

Original comment by robert.c...@gmail.com on 30 Jan 2012 at 10:24

• Changed state: Migrated