← Back to branch summary

~pragmatic-core/pragmatic/scaling_optimisation

« back to all changes in this revision

Viewing changes to python/minimal_example.py

  • Committer: kristian.jensen at ac
  • Date: 2014-01-23 10:35:07 UTC
  • Revision ID: kristian.jensen@ic.ac.uk-20140123103507-xusx40435unfljmt
comments everywhere and advection example

Show diffs side-by-side

added added

removed removed

Lines of Context:
python/minimal_example.py
 
1
### this a testcase for use with DOLFIN/FEniCS and PRAgMaTIc 
 
2
### by Kristian Ejlebjerg Jensen, January 2014, Imperial College London
 
3
### the purpose of the test case is to
 
4
### 1. derive a forcing term that gives a step function solution
 
5
### 2. solve a poisson equation with the forcing term using 2D anisotropic adaptivity
 
6
### 3. calculate and plot the L2error of the resulting solution.
 
7
### the width of the step, the number of solition<->adaptation iterations as well as the
 
8
### error level (eta) are optional input parameters
 
9
 
1
10
from dolfin import *
2
11
from adaptivity2 import metric_pnorm, adapt
3
 
from pylab import hold, show, triplot, tricontourf, colorbar, axis, box, rand, get_cmap
 
12
from pylab import hold, show, triplot, tricontourf, colorbar, axis, box, rand, get_cmap, title
4
13
from pylab import plot as pyplot
5
14
from numpy import array, ones
6
15
import numpy
7
 
set_log_level(WARNING)
8
16
from sympy import Symbol, diff
9
17
from sympy import tanh as pytanh
10
18
from sympy import cos as pysin
11
19
from sympy import sin as pycos
 
20
set_log_level(INFO)
 
21
#parameters["allow_extrapolation"] = True
12
22
 
13
 
def minimal_example(width=5e-2):
 
23
def minimal_example(width=5e-2, Nadapt=10, eta = 0.01):
14
24
    ### CONSTANTS
15
25
    meshsz = 40
16
26
    hd = Constant(width)
17
27
    ### SETUP MESH
18
28
    mesh = RectangleMesh(-0.5,-0.5,0.5,0.5,1*meshsz,1*meshsz,"left/right")
19
 
    ### SETUP SOLUTION
 
29
    ### DERIVE FORCING TERM
20
30
    angle = pi/8 #rand*pi/2
21
 
    #testsol = 'tanh(x[0]/' + str(float(hd)) + ')' #tanh(x[0]/hd)
22
 
    #sx = Symbol('sx'); sy = Symbol('sy'); width_ = Symbol('ww'); angle_ = Symbol('aa')
23
 
    #testsol = pytanh((pycos(angle_)*sx+pysin(angle_)*sy)/width_)    
24
 
    testsol = 'tanh((' + str(cos(angle)) + '*x[0]+'+ str(sin(angle)) + '*x[1])/' + str(float(hd)) + ')' #tanh(x[0]/hd)
25
 
    ddtestsol = str(cos(angle)+sin(angle))+'*2*'+testsol+'*(1-pow('+testsol+',2))/'+str(float(hd)**2)
26
 
    
 
31
    sx = Symbol('sx'); sy = Symbol('sy'); width_ = Symbol('ww'); aa = Symbol('aa')
 
32
    testsol = pytanh((sx*pycos(aa)+sy*pysin(aa))/width_)
 
33
    ddtestsol = str(diff(testsol,sx,sx)+diff(testsol,sy,sy)).replace('sx','x[0]').replace('sy','x[1]')
 
34
    #replace ** with pow
 
35
    ddtestsol = ddtestsol.replace('tanh((x[0]*sin(aa) + x[1]*cos(aa))/ww)**2','pow(tanh((x[0]*sin(aa) + x[1]*cos(aa))/ww),2.)')
 
36
    ddtestsol = ddtestsol.replace('cos(aa)**2','pow(cos(aa),2.)').replace('sin(aa)**2','pow(sin(aa),2.)').replace('ww**2','(ww*ww)')
 
37
    #insert vaulues
 
38
    ddtestsol = ddtestsol.replace('aa',str(angle)).replace('ww',str(width))
 
39
    testsol = str(testsol).replace('sx','x[0]').replace('sy','x[1]').replace('aa',str(angle)).replace('ww',str(width))
 
40
    ddtestsol = "-("+ddtestsol+")"
 
41
    def boundary(x):
 
42
          return x[0]-mesh.coordinates()[:,0].min() < DOLFIN_EPS or mesh.coordinates()[:,0].max()-x[0] < DOLFIN_EPS \
 
43
          or mesh.coordinates()[:,1].min()+0.5 < DOLFIN_EPS or mesh.coordinates()[:,1].max()-x[1] < DOLFIN_EPS  
27
44
    # PERFORM TEN ADAPTATION ITERATIONS
28
 
    for iii in range(10):
 
45
    for iii in range(Nadapt):
29
46
     V = FunctionSpace(mesh, "CG" ,2); dis = TrialFunction(V); dus = TestFunction(V); u = Function(V)
30
 
     R = interpolate(Expression(ddtestsol),V)
31
47
     a = inner(grad(dis), grad(dus))*dx
32
 
     L = R*dus*dx
33
 
     solve(a == L, u, [])
 
48
     L = Expression(ddtestsol)*dus*dx
 
49
     bc = DirichletBC(V, Expression(testsol), boundary)
 
50
     solve(a == L, u, bc)
34
51
     startTime = time()
35
 
     eta = 0.01; H = metric_pnorm(u, mesh, eta, max_edge_ratio=50);   Mp =  project(H,  TensorFunctionSpace(mesh, "CG", 1))
36
 
     mesh = adapt(Mp) 
37
 
     print("total (adapt+metric) time was %0.1fs" % (time()-startTime))
 
52
     H = metric_pnorm(u, eta, max_edge_length=1., max_edge_ratio=50)
 
53
     if iii != Nadapt-1:
 
54
      mesh = adapt(H)
 
55
      L2error = errornorm(Expression(testsol), u, degree_rise=4, norm_type='L2')
 
56
      log(INFO+1,"total (adapt+metric) time was %0.1fs, L2error=%0.0e, nodes: %0.0f" % (time()-startTime,L2error,mesh.num_vertices()))
38
57
    
39
58
    # PLOT MESH
 
59
    testfe = interpolate(u,FunctionSpace(mesh,'CG',1))
40
60
    testf = interpolate(Expression(testsol),FunctionSpace(mesh,'CG',1))
41
61
    vtx2dof = vertex_to_dof_map(FunctionSpace(mesh, "CG" ,1))
42
 
    zz = testf.vector().array()[vtx2dof]; zz[zz==1] -= 1e-16
 
62
    zz = testf.vector().array()[vtx2dof]-testfe.vector().array()[vtx2dof]; zz[zz==1] -= 1e-16
43
63
    hh=tricontourf(mesh.coordinates()[:,0],mesh.coordinates()[:,1],mesh.cells(),zz,100,cmap=get_cmap('binary'))
44
64
    colorbar(hh)
45
65
    
46
66
    hold('on'); triplot(mesh.coordinates()[:,0],mesh.coordinates()[:,1],mesh.cells(),color='r',linewidth=0.5); hold('off')
47
 
    axis('equal'); box('off'); show()
 
67
    axis('equal'); box('off'); title('error'); show()
48
68
 
49
69
if __name__=="__main__":
50
70
 minimal_example()