In [1]:
import os
import sys
import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt
from visualize import drawCoefficient, d3sol, d3solextra
from gridlod import util, world, fem
from gridlod.world import World
import femsolverCoarse, buildcoef2d
For a two dimensional problem, we use the 'buildcoef2d' class. The constructed coefficient can be visualized with a two dimensional plot function that can be found in 'visualize'. For the reference solution we compute a solution on a fine mesh with $1/h=1024$.
In [2]:
bg = 0.01 # background
val = 1 # values
#fine World
fine = 1024
NWorldFine = np.array([fine, fine])
NpFine = np.prod(NWorldFine+1)
NList = [fine,2,4,8,16,32,64,128,256,512]
boundaryConditions = np.array([[0, 0],
[0, 0]])
#coefficient
WormsUltCoefClass = buildcoef2d.Coefficient2d(NWorldFine,
bg = bg,
val = val,
length = 1,
thick = 65,
space = 65,
probfactor = 1,
right = 1,
down = 0,
diagr1 = 0,
diagr2 = 0,
diagl1 = 0,
diagl2 = 0,
LenSwitch = None,
thickSwitch = None,
equidistant = True,
ChannelHorizontal = None,
ChannelVertical = True,
BoundarySpace = None,
Boxes2n = None,
Channels2n = None,
TestExample = True)
A = WormsUltCoefClass.BuildCoefficient()
ABase = A.flatten()
plt.figure("OriginalCoefficient")
drawCoefficient(NWorldFine, ABase)
plt.title("Coefficient")
plt.show()
We compute the FEM approximation for each mesh size and display each solution for comparison. Furthermore, we store the energy error and visualize the error plot it subsequently.
In [3]:
newErrorFine = []
x = []
y = []
for N in NList:
NWorldCoarse = np.array([N,N])
# print "Compute for N = "+ str(N) + " ..."
NCoarseElement = NWorldFine/NWorldCoarse
world = World(NWorldCoarse, NCoarseElement, boundaryConditions)
AFine = fem.assemblePatchMatrix(NWorldFine, world.ALocFine, ABase)
#grid nodes
xpCoarse = util.pCoordinates(NWorldCoarse).flatten()
NpCoarse = np.prod(NWorldCoarse+1)
f = np.ones(NpCoarse)
uCoarseFull = femsolverCoarse.solveCoarse_fem(world, ABase, f, boundaryConditions)
basis = fem.assembleProlongationMatrix(NWorldCoarse, NCoarseElement)
uLodCoarse = basis*uCoarseFull
if N == fine:
uSol = uCoarseFull
else:
newErrorFine.append(np.sqrt(np.dot(uSol - uLodCoarse, AFine*(uSol - uLodCoarse))))
x.append(N)
y.append(1./N)
ymin = np.min(uSol)
ymax = np.max(uSol)
if np.size(x)==2:
fig = plt.figure('FEM-Solution')
fig.subplots_adjust(left=0.01,bottom=0.04,right=0.99,top=0.95,wspace=0,hspace=0.2)
ax = fig.add_subplot(241, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
elif np.size(x)==3:
ax = fig.add_subplot(242, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
elif np.size(x)==4:
ax = fig.add_subplot(243, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
elif np.size(x)==5:
ax = fig.add_subplot(244, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
elif np.size(x)==6:
ax = fig.add_subplot(245, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
elif np.size(x)==7:
ax = fig.add_subplot(246, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
elif np.size(x)==8:
ax = fig.add_subplot(247, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
elif np.size(x)==9:
ax = fig.add_subplot(248, projection='3d')
ax.set_title('1/h= ' + str(N),fontsize="small")
if N == fine:
d3sol(NWorldCoarse, uCoarseFull)
elif N == 2:
init = 0
else:
d3solextra(NWorldCoarse, uCoarseFull, fig, ax, ymin, ymax)
plt.show()
In [4]:
plt.figure("Error")
plt.loglog(x,newErrorFine,'o-', basex=2, basey=2)
plt.loglog(x,y,'--k',basex=2, basey=2, linewidth=1, alpha=0.3)
plt.ylabel('Energy error')
plt.xlabel('$1/h$')
plt.subplots_adjust(left=0.1,bottom=0.1,right=0.95,top=0.95,wspace=0.2,hspace=0.2)
plt.title('Energy error for the standard FEM')
plt.grid(True,which="both",ls="--")
plt.show()