diff --git a/1d/main.py b/1d/main.py
new file mode 100644
index 0000000000000000000000000000000000000000..5644f258ee108d0fe9dde79c0dc43672c5012413
--- /dev/null
+++ b/1d/main.py
@@ -0,0 +1,196 @@
+import myfun as my
+import numpy as np
+from matplotlib import pyplot as plt
+
+# CHOOSE METHOD
+method = 'wasserstein' #'wasserstein' or 'upwinding'
+
+T = 30*10**(-3) # time interval [0,T]
+
+Nx = 103 # choose fineness of spatial discretization
+Nt = 200 # choose fineness of temporal discretization
+hx = 1/(Nx+2) # assumption equidistant mesh, spatial mesh size
+ht = (T)/(Nt+1) # assumption equidistant steps, temporal mesh size
+
+print('hx',hx)
+print('ht',ht)
+
+faces = np.linspace(0,1,Nx+3) # Cell interfaces x_(i+1/2) for i=0,1,...,N and boundary faces x_-1/2 = 0, x_N+1/2 = 1.
+# +3 as we have N+1 interfaces and 2 boundary faces
+faces = np.array(faces) # [-1, 0,1,2,...,N ,N+1] -> length = Nx+3
+
+# compute array of midpoints
+midpoints = []
+for i in range(Nx+2) :
+ midpoints.append((faces[i]+faces[i+1])/2) # cell midpoints
+midpoints = np.array(midpoints) # -> length = Nx+2
+
+# pre-allocate
+rho = np.zeros([Nt,Nx+2])
+cc = np.zeros([Nt,Nx+2])
+
+# CHOOSE INITIAL CONDITION
+rho_0 = lambda x: 1.3*np.sin(np.pi*x)*np.exp(-25*(x-1/2)**2) # pre-factors scales measure to sum up to 1 over midpoints
+# rho_0 = lambda x: (np.sin(np.pi*x)*np.exp(-50*(x-1/4)**2) + np.sin(np.pi*x)*np.exp(-50*(x-3/4)**2)) # double gaussian
+
+# discrete initial condition via evalutation at midpoints, i.e. projection to discrete space
+rho_h0 = rho_0(midpoints) # initial value t = 0, N+3 faces result in N+2 elements, need to sum to 1 (represents measure), may not start with uniform distribution
+
+# save initial discrete rho0
+rho[0][:] = rho_h0 # t=0
+
+# initialize object
+discr_E = []
+
+
+if method == 'upwinding' :
+
+ for n in range(Nt-1) : # go through time intervals
+
+ # GET C FROM RHO
+ c = my.getc(rho[n,:],midpoints,faces,hx) # paramaters to get function c_h^n(x) via c_h^n(x) = sum_i (c_i*hat_i(x))
+ cc[n,:] = c
+
+ # get discrete energy
+ discr_E.append(my.discrete_energy(c,rho[n,:],hx))
+ rhs = []
+
+ # APPLY NUMERICAL SCHEME
+ for i in range(Nx+2) : # assumption gamma = 1, len(rho) = Nx+2, going over midpoints
+ #print(i)
+ if i == Nx+2-1 : # periodic boundary data
+ dxch_irechts = c[i]*my.dxhat(i,midpoints,hx,faces[0])+c[0]*my.dxhat(0,midpoints,hx,faces[0])
+ dxch_ilinks = c[i-1]*my.dxhat(i-1,midpoints,hx,faces[i])+c[i]*my.dxhat(i,midpoints,hx,faces[i]) #x_{i-1/2}
+
+ if dxch_irechts >= 0 :
+ if dxch_ilinks >= 0 :
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][i]-dxch_ilinks*rho[n][i-1]))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][i]-dxch_ilinks*rho[n][i]))
+ else :
+ if dxch_ilinks >= 0:
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][0]-dxch_ilinks*rho[n][i-1]))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][0]-dxch_ilinks*rho[n][i]))
+
+ else : # interior faces
+ dxch_irechts = c[i]*my.dxhat(i,midpoints,hx,faces[i+1])+c[i+1]*my.dxhat(i+1,midpoints,hx,faces[i+1])
+ dxch_ilinks = c[i-1]*my.dxhat(i-1,midpoints,hx,faces[i])+c[i]*my.dxhat(i,midpoints,hx,faces[i])
+
+ if dxch_irechts >= 0 :
+ if dxch_ilinks >= 0:
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][i]-dxch_ilinks*rho[n][i-1]))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][i]-dxch_ilinks*rho[n][i]))
+ else :
+ if dxch_ilinks >= 0:
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][i+1]-dxch_ilinks*rho[n][i-1]))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_irechts*rho[n][i+1]-dxch_ilinks*rho[n][i]))
+
+ rho[n+1][:] = my.update_in_time(rhs,Nx,ht,hx)
+
+ discr_E = np.array(discr_E)
+
+
+elif method == 'wasserstein' :
+
+ for n in range(Nt-1) : # go through time intervals
+
+ # GET C FROM RHO
+ c = my.getc(rho[n,:],midpoints,faces,hx) # paramaters to get function c_h^n(x) via c_h^n(x) = sum_i (c_i*hat_i(x))
+ cc[n,:] = c
+
+ # get discrete energy
+ discr_E.append(my.discrete_energy(c,rho[n,:],hx))
+ rhs = []
+
+ # APPLY NUMERICAL SCHEME
+ for i in range(Nx+2) : # assumption gamma = 1, len(rho) = Nx+2
+ if i == Nx+2-1 : # periodic boundary data
+ dxch_iplus = c[i]*my.dxhat(i,midpoints,hx,faces[i+1])+c[0]*my.dxhat(0,midpoints,hx,faces[i+1])
+ dxch_iminus = c[i-1]*my.dxhat(i-1,midpoints,hx,faces[i])+c[i]*my.dxhat(i,midpoints,hx,faces[i])
+ if np.abs(rho[n][i]-rho[n][i-1]) < 10**-6 or rho[n][i] < 10**-6 or rho[n][i-1] < 10**-6 :
+ if np.abs(rho[n][0]-rho[n][i]) < 10**-6 or rho[n][i] < 10**-6 or rho[n][0] < 10**-6 :
+ rhs.append(rho[n][i] - ht/hx*(dxch_iplus*1/2*(rho[n][0]+rho[n][i])-dxch_iminus*1/2*(rho[n][i]+rho[n][i-1])))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_iplus*(rho[n][0]-rho[n][i])/(np.log(rho[n][0])-np.log(rho[n][i]))-dxch_iminus*1/2*(rho[n][i]+rho[n][i-1])))
+ else :
+ if np.abs(rho[n][0]-rho[n][i]) < 10**-6 or rho[n][i] < 10**-6 or rho[n][0] < 10**-6 :
+ rhs.append(rho[n][i] - ht/hx*1/2*(dxch_iplus*(rho[n][0]+rho[n][i])-dxch_iminus*(rho[n][i]-rho[n][i-1])/(np.log(rho[n][i])-np.log(rho[n][i-1]))))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_iplus*(rho[n][0]-rho[n][i])/(np.log(rho[n][0])-np.log(rho[n][i]))-dxch_iminus*(rho[n][i]-rho[n][i-1])/(np.log(rho[n][i])-np.log(rho[n][i-1]))))
+
+ else : # interior faces
+ dxch_iplus = c[i]*my.dxhat(i,midpoints,hx,faces[i+1])+c[i+1]*my.dxhat(i+1,midpoints,hx,faces[i+1])
+ dxch_iminus = c[i-1]*my.dxhat(i-1,midpoints,hx,faces[i])+c[i]*my.dxhat(i,midpoints,hx,faces[i])
+ if np.abs(rho[n][i]-rho[n][i-1]) < 10**-6 or rho[n][i] < 10**-6 or rho[n][i-1] < 10**-6 :
+ if np.abs(rho[n][i+1]-rho[n][i]) < 10**-6 or rho[n][i] < 10**-6 or rho[n][i+1] < 10**-6 :
+ rhs.append(rho[n][i] - ht/hx*(dxch_iplus*1/2*(rho[n][i+1]+rho[n][i])-dxch_iminus*1/2*(rho[n][i]+rho[n][i-1])))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_iplus*(rho[n][i+1]-rho[n][i])/(np.log(rho[n][i+1])-np.log(rho[n][i]))-dxch_iminus*1/2*(rho[n][i]+rho[n][i-1])))
+ else :
+ if np.abs(rho[n][i+1]-rho[n][i]) < 10**-6 or rho[n][i] < 10**-6 or rho[n][i+1] < 10**-6 :
+ rhs.append(rho[n][i] - ht/hx*1/2*(dxch_iplus*(rho[n][i+1]+rho[n][i])-dxch_iminus*(rho[n][i]-rho[n][i-1])/(np.log(rho[n][i])-np.log(rho[n][i-1]))))
+ else :
+ rhs.append(rho[n][i] - ht/hx*(dxch_iplus*(rho[n][i+1]-rho[n][i])/(np.log(rho[n][i+1])-np.log(rho[n][i]))-dxch_iminus*(rho[n][i]-rho[n][i-1])/(np.log(rho[n][i])-np.log(rho[n][i-1]))))
+
+ rho[n+1][:] = my.update_in_time(rhs,Nx,ht,hx)
+
+ # PLOT FUNCTION PROFILES AS SNAPSHOTS AT TIME STEP n
+ # does not work for all Nx
+ if n == 42:
+ c_h = lambda x : my.cfunc(cc[n][:],midpoints,hx,x)
+ lin_rho = lambda x : my.linear_interpol(rho[n][:],midpoints,hx,x)
+ const_rho = lambda x: my.const_rho(rho[n][:],midpoints,hx,x)
+
+ #compute function values
+ xi = np.linspace(0,1,300)
+ val1 = []
+ val2 = []
+ val3 = []
+ for x in xi:
+ val1.append(const_rho(x))
+ val2.append(lin_rho(x))
+ val3.append(c_h(x))
+
+ plt.plot(xi,val1,'steelblue',xi,val2,'deepskyblue')
+ plt.plot(midpoints,np.zeros(len(midpoints)),'r.')
+ plt.xlabel('x')
+ plt.ylabel('bacterial density')
+ plt.legend(['finite volume sol.','lin. interpolation','midpoints x_i'])
+ time = n*ht
+ title = 'bacterial density at time '+'%f' % time
+ plt.title(title)
+ plt.show()
+
+ plt.plot(xi,val3,'green')
+ plt.plot(midpoints,np.zeros(len(midpoints)),'r.')
+ plt.xlabel('x')
+ plt.ylabel('chemical concentration')
+ title = 'chemical concentration at time '+'%f' % time
+ plt.title(title)
+ plt.legend(['finite element sol.','midpoints x_i'])
+ plt.show()
+
+
+ # save discrete energy for this time step
+ discr_E = np.array(discr_E)
+
+else :
+ print('Error: Wrong string input')
+
+
+# HEATMAP PLOT OF EVOLUTION
+plt.pcolormesh(rho,cmap = 'jet')
+plt.xlabel('space')
+plt.ylabel('time')
+plt.colorbar()
+plt.show()
+
+# PLOT DISCRETE ENERGY
+plt.plot(np.linspace(0,T,Nt-1),discr_E)
+plt.xlabel('time')
+plt.ylabel('discrete Energy')
+plt.title('Discrete energy dissipation')
+plt.show()
diff --git a/1d/myfun.py b/1d/myfun.py
new file mode 100644
index 0000000000000000000000000000000000000000..ed371dfcedb9f53fae09a09da7db2fe92082d2e7
--- /dev/null
+++ b/1d/myfun.py
@@ -0,0 +1,315 @@
+import numpy as np
+from scipy.sparse import csr_matrix
+from scipy.sparse.linalg import spsolve
+
+
+# FINITE ELEMENT
+
+def linear_interpol(rho,midpoints,h,x) :
+ # Input : current approximation rho^n, midpoints, mesh size h, point x
+ # output : value of linear interpolation at x
+
+ # find element in which x lies use midpoint[i] = h(i+1/2)
+
+ i0 = int(np.floor((x/h-1/2)))
+ i1 = int(np.ceil((x/h-1/2)))
+
+ if i0 == -1 : # periodic boundary condition
+ y1 = rho[i1]
+ x1 = midpoints[i1]
+ y0 = 1/2*(rho[0]+rho[-1])
+ x0 = 0
+
+ #get linear interpolation
+ m = (y1-y0)/(x1-x0) # /h
+ b = y1-m*x1
+
+ y = m*x+b
+
+ elif i1 == len(midpoints) : # periodic boundary condition
+ y1 = 1/2*(rho[0]+rho[-1])
+ x1 = 1
+ y0 = rho[i0]
+ x0 = midpoints[i0]
+
+ #get linear interpolation
+ m = (y1-y0)/(x1-x0) # /h
+ b = y1-m*x1
+
+ y = m*x+b
+
+ else :
+
+ if i0 == i1 : # falls rundungsfehler
+ i0 -= 1
+
+ y1 = rho[i1]
+ x1 = midpoints[i1]
+ y0 = rho[i0]
+ x0 = midpoints[i0]
+
+ #get linear interpolation
+ m = (y1-y0)/(x1-x0) # /h
+ b = y1-m*x1
+
+ y = m*x+b
+
+ return y
+
+
+def hat(i,midpoints,h,x) :
+
+ values = np.zeros(len(midpoints))
+ values[i] = 1
+
+ y = linear_interpol(values,midpoints,h,x)
+
+ return y
+
+
+def dxhat(i,midpoints,h,x) :
+
+ values = np.zeros(len(midpoints))
+ values[i] = 1
+
+ # find element in which x lies use midpoint[i] = h(i+1/2)
+ i0 = int(np.floor((x/h-1/2)))
+ i1 = int(np.ceil((x/h-1/2)))
+
+ if i0 == -1: # periodic boundary condition
+ y1 = values[i1]
+ x1 = midpoints[i1]
+ y0 = 1/2*(values[0]+values[-1])
+ x0 = 0
+
+ #get linear interpolation
+ m = (y1-y0)/(x1-x0) # /h
+
+
+ elif i1 == len(midpoints) : # periodic boundary condition
+ y1 = 1/2*(values[0]+values[-1])
+ x1 = 1
+ y0 = values[i0]
+ x0 = midpoints[i0]
+
+ #get linear interpolation
+ m = (y1-y0)/(x1-x0) # /h
+
+ else :
+
+ y1 = values[i1]
+ x1 = midpoints[i1]
+ y0 = values[i0]
+ x0 = midpoints[i0]
+
+ #get linear interpolation
+ m = (y1-y0)/(x1-x0) # /h
+
+
+ return m
+
+
+def quadrature(I,g) :
+ # inupt: interval I = [a,b]
+
+ areaI = I[1]-I[0]
+
+ trafo = lambda y : y*(I[1]-I[0])/2+(I[1]+I[0])/2 # y in [-1,1]
+
+ xi = np.array([-np.sqrt(1/3),np.sqrt(1/3)]) #on [-1,1]
+ #xi = [-1,1]
+ wi = [1,1]
+
+ # xi = [0]
+ # wi = [2]
+
+ # xi = np.array([-np.sqrt(3/5),0,np.sqrt(3/5)]) #on [-1,1]
+ # wi = [5/9,8/9,5/9]
+
+ val = 0
+ for k in range(len(wi)) :
+ val += wi[k]*g(trafo(xi[k])) # get integral
+
+ avrg = val/2 # get average
+ integral = avrg*areaI
+
+ return integral
+
+def midpoint_rule(I,g) :
+ # inupt: interval I = [a,b]
+
+ areaI = I[1]-I[0]
+
+ trafo = lambda y : y*(I[1]-I[0])/2+(I[1]+I[0])/2 # y in [-1,1]
+
+ xi = [0]
+ wi = [2]
+
+ val = 0
+ for k in range(len(wi)) :
+ val += wi[k]*g(trafo(xi[k])) # get integral
+
+ avrg = val/2 # get average
+ integral = avrg*areaI
+
+ return integral
+
+def trapezoidal_rule(I,g) :
+ # inupt: interval I = [a,b]
+
+ areaI = I[1]-I[0]
+
+ trafo = lambda y : y*(I[1]-I[0])/2+(I[1]+I[0])/2 # y in [-1,1]
+
+ xi = [-1,1]
+ wi = [1,1]
+
+ val = 0
+ for k in range(len(wi)) :
+ val += wi[k]*g(trafo(xi[k])) # get integral
+
+ avrg = val/2 # get average
+ integral = avrg*areaI
+
+ return integral
+
+
+def getc(rho,midpoints,faces,h) :
+
+
+ lin_rho = lambda x : linear_interpol(rho,midpoints,h,x)
+
+ row = []
+ col = []
+ data1 = []
+ data2 = []
+ b = []
+
+ for i in range(len(midpoints)) :
+
+ f = lambda x: hat(i,midpoints,h,x)*lin_rho(x)
+
+ # if i == (len(midpoints)-1):
+ # rhs = quadrature([midpoints[i-1],1],f) + quadrature([0,midpoints[1]],f) # periodic boundary
+ # elif i == 0 :
+ # rhs = quadrature([0,midpoints[i+1]],f) + quadrature([midpoints[-1],1],f) # periodic boundary
+ # else :
+ # rhs = quadrature([midpoints[i-1],midpoints[i+1]],f)
+
+ rhs = midpoint_rule([faces[i],faces[i+1]],f)
+
+ for j in range(len(midpoints)) :
+ if abs(i-j) < 2 :
+ g = lambda x: hat(i,midpoints,h,x)*hat(j,midpoints,h,x)
+ dg = lambda x: dxhat(i,midpoints,h,x)*dxhat(j,midpoints,h,x)
+
+ if i == (len(midpoints)-1) :
+ val2 = trapezoidal_rule([midpoints[i],1],g) + trapezoidal_rule([0,midpoints[0]],g)
+ # val2 = midpoint_rule([midpoints[i],1],g) + midpoint_rule([0,midpoints[0]],g)
+
+ else :
+
+ val2 = trapezoidal_rule([midpoints[i],midpoints[i+1]],g) # for A2
+ # val2 = midpoint_rule([midpoints[i],midpoints[i+1]],g) # for A1
+
+ val1 = trapezoidal_rule([faces[i],faces[i+1]],dg)
+ # val2 = trapezoidal_rule([midpoints[i],midpoints[i+1]],g)
+
+ row.append(i)
+ col.append(j)
+ data1.append(val1)
+ data2.append(val2)
+ b.append(rhs)
+
+
+
+ sparseA1 = csr_matrix((data1, (row, col)), shape = (len(midpoints), len(midpoints)))#.toarray()
+ sparseA2 = csr_matrix((data2, (row, col)), shape = (len(midpoints), len(midpoints)))#.toarray()
+
+ sparseA = sparseA1 + sparseA2
+
+ # plt.spy(sparseA2.toarray())
+ # plt.show()
+
+ c = spsolve(sparseA, b, permc_spec=None, use_umfpack=True) # solve system of equations
+ # c_h^n(x) = sum_i (c_i*hat_i(x))
+ return c
+
+
+def cfunc(cc,midpoints,h,x) :
+
+ i0 = int(np.floor((x/h-1/2)))
+ i1 = int(np.ceil((x/h-1/2)))
+
+ if i1 == len(midpoints) :
+ i1 = 0
+
+ # print(ii)
+ val = 0
+ for i in [i0,i1] :
+ val += cc[i]*hat(i,midpoints,h,x)
+
+ return val
+
+
+def const_rho(rho,midpoints,h,x) :
+
+ i = int(np.round((x/h-1/2)))
+ val = rho[i]
+
+ return val
+
+
+def update_in_time(rhs,Nx,ht,hx) :
+
+ row = []
+ col = []
+ data = []
+
+ row.append(0) # periodic boundary condition
+ col.append(Nx+1)
+ data.append(-ht/(hx**2))
+
+ for i in range(Nx+2) :
+ for j in range(Nx+2) :
+ if abs(i-j)<2:
+ if i == j :
+ val = 1+(2*ht)/(hx**2)
+ else :
+ val = -ht/(hx**2)
+
+ row.append(i)
+ col.append(j)
+ data.append(val)
+
+ row.append(Nx+1) # periodic boundary condition
+ col.append(0)
+ data.append(-ht/(hx**2))
+
+ sparseAt = csr_matrix((data, (row, col)), shape = (Nx+2, Nx+2))#.toarray()
+ # plt.spy(csr_matrix((data, (row, col)), shape = (Nx+2, Nx+2)).toarray())
+ # plt.show()
+ rho_time_update = spsolve(sparseAt, rhs, permc_spec=None, use_umfpack=True) # solve system of equations
+
+ #print(len(rho_time_update))
+ return rho_time_update
+
+
+def discrete_energy(c,rho,hx) :
+
+ # input: discrete values rho and c, and spatial step size hx
+ # discrete energey w.r.t. rho and c
+
+ value = 0
+
+ for i in range(len(rho)) :
+
+ # MIDPOINT RULE
+ value += hx*(rho[i]*np.log(rho[i])-1/2*rho[i]*c[i]) # cc[i] is fine as hat(i) = 1 in midpoint[i]
+
+
+ return value
+
+
+
+
diff --git a/2d/__pycache__/myfun.cpython-312.pyc b/2d/__pycache__/myfun.cpython-312.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..39b78f767f2996c22ecb00cf1bd39265a7a5aea7
Binary files /dev/null and b/2d/__pycache__/myfun.cpython-312.pyc differ
diff --git a/2d/compute_errorestimates.py b/2d/compute_errorestimates.py
new file mode 100644
index 0000000000000000000000000000000000000000..c4e7adefd1f2dd35fba50fb6a041371a1280b52f
--- /dev/null
+++ b/2d/compute_errorestimates.py
@@ -0,0 +1,252 @@
+import myfun as my
+import numpy as np
+import time
+import pickle
+
+tic = time.time()
+
+# ---------------------------------------------------- COMPUTE A POSTERIORI ERROR ESTIMATES -----------------------------------------------------
+# This script computes the individual terms of the a posteriori error estimates \theta_K for each element K \in \mathca{T}_h
+
+# SETTINGS FOR SCHEME
+nodal_avrg_choice = 'arithmetic-mean'
+method = 'wasserstein' # using log-mean convective flux
+boundary = 'Neumann'
+interpol = 'morley'
+
+test = 'manufactured1' # 'diff', 'blow-up', 'manufactured' or 'manufactured1'
+fineness = 0 # number of refinements of mesh before calculating numerical solution
+Nt = 8 # number of time steps
+
+
+# initialize different test cases
+if test == 'blow-up' :
+# CONVECTION DEOMINATED REGIME i.e. BLOW-UP
+
+ def rho0(x) :
+ val = 10**(3)*np.exp(-((x[0]-0.5)**2+(x[1]-0.5)**2)/10**(-2))
+ return val
+
+ T = 0.0045 # time interval [0,T]
+
+elif test == 'diff' :
+# DIFFUSION DOMINATED REGIME
+ def rho0(x) :
+ val = 1.3*np.exp(-25*(x[0]-1/2)**2-25*(x[1]-1/2)**2)
+ return val
+
+ T = 0.05 # time interval [0,T]#
+
+elif test == 'manufactured' :
+# MANUFACTURED SOLUTION WITH T = 0.5
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ T = 0.5 # time interval [0,T]#
+
+elif test == 'manufactured1' :
+# MANUFACTURED SOLUTION WITH T = 3
+
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ T = 3 # time interval [0,T]#
+
+else :
+ print('Warning: Wrong string input for test.')
+
+
+# SCHEME/PROBLEM PARAMETERS
+hx = 1/2
+ht = (T)/(Nt)
+
+
+# initialize zerost mesh for FV approximations
+[tri,points,nppoints,numtri] = my.initialmesh()
+K = nppoints[tri.simplices]
+
+for i in range(fineness) :
+ hx = 1/2*hx
+ [tri,points,nppoints,numtri] = my.refinemesh(points, K, numtri)
+ K = nppoints[tri.simplices]
+
+
+#get dual mesh for FE method to approximate chemical density
+[K_dual, tri_dual, points_dual, nppoints_dual, numtri_dual] = my.getdualmesh(points, tri, fineness, K)
+K_dual = np.array(K_dual)
+[tri_dual,K_dual] = my.orientation_dualmesh(tri_dual,K_dual) # enforce positive orientation of all triangles in dual mesh
+
+# compute mesh topology for faster find_simplex routine for dual mesh
+[oriented_edges,structure,edge_to_index] = my.get_edges_structure(K_dual)
+my.init_indexE(numtri_dual)
+
+print('Compute error estimates')
+print(interpol+'_'+boundary+'_'+method+'_'+test)
+print('fineness',fineness)
+print('numtri',numtri)
+print('hx',hx)
+print('ht',ht)
+print('numtri_dual',numtri_dual)
+print('Nt',Nt)
+
+max1 = 0 # need this for more efficiency calculating theta_omega
+max2 = 0 # need this for more efficiency calculating theta_omega
+for n in range(Nt) : # range(Nt) :
+
+ progress = n/(Nt)*100
+ print("%.2f" % progress, 'procent of progress made.')
+
+ if n == 0: # slightly different setting at early times
+
+ #GET ARTIFICIAL MORLEY0
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_c at time step'+str(n)+'.p'
+ cc_0 = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # compute c related objects
+ v_0 = lambda x: my.grad_cfunc(cc_0,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) # piecewiese constant
+ [lin_vx_0,lin_vy_0] = my.get_linv(K_dual,nppoints_dual,numtri_dual,tri_dual,v_0,nodal_avrg_choice) # get linear coefficients
+ Laplace_c_0 = lambda x : my.approx_Laplace_cn(lin_vx_0,lin_vy_0,x,oriented_edges,structure,edge_to_index) # approximation of Laplace(c)=div(v)
+
+ v_m1 = lambda x : v_0(x) # as morley0 uses c0 instaed of c^(n-1) in n = 0 case
+
+ # load data
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n+1)+'.p'
+ [ht_0,rho_p1,qq0_p1,betaKE_p1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n)+'.p'
+ [rho_0,rho0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # get artificial morley interpolation
+ [qq0_0,betaKE_0] = my.get_morley0_coeff(rho_0,nppoints,numtri,tri,K)
+
+ #dummy variables
+ rho_m1 = 0*rho_0
+ ht_m1 = 0
+
+ # COMPUTE A POSTERIORI ERROR ESTIMATES
+ theta_early = my.get_early_time_error(K,tri,v_0,qq0_0,qq0_p1,betaKE_0,betaKE_p1)
+ print('extra term for early times computed.')
+ theta_res0 = my.get_theta_res(cc_0,Laplace_c_0,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index)
+ print('elementwise residual computed.')
+ theta_diff0 = my.get_theta_diff(qq0_0,betaKE_0,K,tri)
+ print('grad morley jump terms computed.')
+ [theta_time0,theta_time1] = my.get_theta_time(qq0_0,qq0_p1,betaKE_0,betaKE_p1,rho_m1,rho_0,rho_p1,ht_m1,ht_0,K,tri,n)
+ print('time contributions computed.')
+ [theta_Omega,max1,max2] = my.get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,n)
+ print('global convective term computed.')
+
+ # SAVE A POSTERIORI ERROR ESTIMATES
+ data = []
+ data.append(theta_early)
+ data.append(theta_res0)
+ data.append(theta_diff0)
+ data.append(theta_time0)
+ data.append(theta_Omega)
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ pickle.dump(data,open('pickle_files/'+pickle_name,'wb')) # store data
+
+
+ else :
+
+ if n == 1:
+
+ # time-update objects
+ cc_m1 = cc_0
+ v_m1 = v_0
+ Laplace_c_m1 = Laplace_c_0
+ rho_m1 = rho_0
+ rho_0 = rho_p1
+ qq0_0 = qq0_p1
+ betaKE_0 = betaKE_p1
+ ht_m1 = ht_0
+
+ # load data
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_c at time step'+str(n)+'.p'
+ cc_0 = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n+1)+'.p'
+ [ht_0,rho_p1,qq0_p1,betaKE_p1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # compute c related objects
+ v_0 = lambda x: my.grad_cfunc(cc_0,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) # piecewiese constant
+ [lin_vx_0,lin_vy_0] = my.get_linv(K_dual,nppoints_dual,numtri_dual,tri_dual,v_0,nodal_avrg_choice) # get linear coefficients
+ Laplace_c_0 = lambda x : my.approx_Laplace_cn(lin_vx_0,lin_vy_0,x,oriented_edges,structure,edge_to_index)
+
+ # dummy variables
+ betaKE_m1 = 0*np.array(betaKE_0)
+ v_m2 = lambda x : 0*np.array(v_m1(x))
+
+
+ else :
+
+ # time-update objects
+ v_m2 = v_m1
+ v_m1 = v_0
+ Laplace_c_m1 = Laplace_c_0
+ rho_m1 = rho_0
+ rho_0 = rho_p1
+ qq0_0 = qq0_p1
+ betaKE_m1 = betaKE_0
+ betaKE_0 = betaKE_p1
+ ht_m1 = ht_0
+
+ # load data
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_c at time step'+str(n)+'.p'
+ cc_0 = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n+1)+'.p'
+ [ht_0,rho_p1,qq0_p1,betaKE_p1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # compute c related objects
+ v_0 = lambda x: my.grad_cfunc(cc_0,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) # piecewiese constant
+ [lin_vx_0,lin_vy_0] = my.get_linv(K_dual,nppoints_dual,numtri_dual,tri_dual,v_0,nodal_avrg_choice) # get linear coefficients
+ Laplace_c_0 = lambda x : my.approx_Laplace_cn(lin_vx_0,lin_vy_0,x,oriented_edges,structure,edge_to_index)
+
+
+ # COMPUTE A POSTERIORI ERROR ESTIMATES
+ theta_res0 = my.get_theta_res(cc_m1,Laplace_c_m1,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index)
+ print('elementwise residual computed.')
+ theta_diff0 = my.get_theta_diff(qq0_0,betaKE_0,K,tri)
+ print('grad morley jump terms computed.')
+ [theta_time0,theta_time1] = my.get_theta_time(qq0_0,qq0_p1,betaKE_0,betaKE_p1,rho_m1,rho_0,rho_p1,ht_m1,ht_0,K,tri,n)
+ print('time contributions computed.')
+ conv_terms = my.get_conv_terms(v_0,qq0_p1,betaKE_p1,rho_p1,K,tri)
+ print('convective term computed.')
+ [theta_Omega,max1,max2] = my.get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,n)
+ print('global convective term computed.')
+
+ # SAVE A POSTERIORI ERROR ESTIMATES
+ data = []
+ data.append(theta_res0)
+ data.append(theta_diff0)
+ data.append(theta_time0)
+ data.append(theta_time1)
+ data.append(conv_terms)
+ data.append(theta_Omega)
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ pickle.dump(data,open('pickle_files/'+pickle_name,'wb')) # store data
+
+progress = 100
+print("%.2f" % progress, 'procent of progress made.')
+elapsed = time.time() - tic
+print('This took',"%.2f" % round(elapsed/60, 2), 'minutes.')
+
+
diff --git a/2d/exact_error.py b/2d/exact_error.py
new file mode 100644
index 0000000000000000000000000000000000000000..95d8e57264f7e148daee277c878e45bdab389b4f
--- /dev/null
+++ b/2d/exact_error.py
@@ -0,0 +1,201 @@
+import myfun as my
+import numpy as np
+import time
+import pickle
+
+tic = time.time()
+
+# ----------------------------------------------- COMPUTE "EXACT" ERROR OF MANUFACTURED SOLUTION -----------------------------------------------------
+# This script computes the LinfL2-L2H1 norm of the "exact" error subject to manufactured solutions.
+
+# SETTINGS FOR SCHEME
+nodal_avrg_choice = 'least-squares'
+method = 'wasserstein' # using log-mean convective flux
+boundary = 'Neumann'
+interpol = 'morley'
+
+test = 'manufactured1' # 'diff' or 'blow-up' or 'manufactured'
+fineness = 3 # number of refinements of mesh before calculating numerical solution
+Nt = 8 # number of time steps, [1, 2, 3, 4, 6 ,8 , 16, 32, 64]
+
+# initialize different test cases; only works for manufactured solutions
+if test == 'manufactured' :
+# MANUFACTURED SOLUTION WITH T = 0.5
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def grad_exact_rho(x,t) :
+ val = -1.3*120/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*np.array([x[0]-1/2,x[1]-1/2])
+ return val
+
+ T = 0.5 # time interval [0,T]#
+
+elif test == 'manufactured1' :
+# MANUFACTURED SOLUTION WITH T = 3
+
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def grad_exact_rho(x,t) :
+ val = -1.3*120/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*np.array([x[0]-1/2,x[1]-1/2])
+ return val
+
+ T = 3 # time interval [0,T]#
+
+else :
+ print('Warning: Wrong string input for test.')
+
+
+# SCHEME/PROBLEM PARAMETERS
+hx = 1/2
+ht = (T)/(Nt)
+
+# initialize zerost mesh for FV approximations
+[tri,points,nppoints,numtri] = my.initialmesh()
+K = nppoints[tri.simplices]
+
+# refine zerost mesh subject to fineness
+for i in range(fineness) :
+ hx = 1/2*hx
+ [tri,points,nppoints,numtri] = my.refinemesh(points, K, numtri)
+ K = nppoints[tri.simplices]
+
+
+#get dual mesh for FE method to approximate chemical density
+[K_dual, tri_dual, points_dual, nppoints_dual, numtri_dual] = my.getdualmesh(points, tri, fineness, K)
+K_dual = np.array(K_dual)
+[tri_dual,K_dual] = my.orientation_dualmesh(tri_dual,K_dual) # enforce positive orientation of all triangles in dual mesh
+
+# compute mesh topology for faster find_simplex routine for dual mesh
+[oriented_edges,structure,edge_to_index] = my.get_edges_structure(np.array(K_dual))
+my.init_indexE(numtri_dual)
+
+print('Compute exact error')
+print(interpol+'_'+boundary+'_'+method+'_'+test)
+print('fineness',fineness)
+print('numtri',numtri)
+print('hx',hx)
+print('ht',ht)
+print('numtri_dual',numtri_dual)
+print('Nt',Nt)
+
+# inititlize objects
+L2_array = []
+L2H1_sum = 0
+
+for n in range(Nt) : # range(Nt) :
+
+ progress = n/(Nt)*100
+ print("%.2f" % progress, 'procent of progress made.')
+
+ if n == 0 : # slightly different setting at early times
+
+ # load data
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n)+'.p'
+ [rho_0,rho0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n+1)+'.p'
+ [ht_0,rho_p1,qq0_p1,betaKE_p1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ #GET ARTIFICIAL MORLEY0
+ [qq0_0,betaKE_0] = my.get_morley0_coeff(rho_0,nppoints,numtri,tri,K)
+
+ else :
+
+ if n == 1:
+
+ # time-update objects
+ qq0_0 = qq0_p1
+ betaKE_0 = betaKE_p1
+
+ # load data
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n+1)+'.p'
+ [ht_0,rho_p1,qq0_p1,betaKE_p1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+
+ else :
+
+ # time-update objects
+ qq0_0 = qq0_p1
+ betaKE_0 = betaKE_p1
+
+ # load data
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n+1)+'.p'
+ [ht_0,rho_p1,qq0_p1,betaKE_p1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+
+ # COMPUTE RECONSTRUCTION OF MORLEY-TYPE FOR BOTH TIME STEPS, INCLUDUNG ITS GRAD
+ morley_interpol0 = lambda x : my.morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=-1) # use rho0 here?
+ grad_morley_interpol0 = lambda x : my.grad_morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=-1)
+ morley_interpol1 = lambda x : my.morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=-1)
+ grad_morley_interpol1 = lambda x : my.grad_morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=-1)
+
+ # COMPUTE FULL RECONSTRUCTION (time and space), INCLUDING ITS GRAD
+ morley_spacetime = lambda x,t : morley_interpol1(x)*(t-n*ht)/((n+1)*ht-n*ht) - morley_interpol0(x)*(t-(n+1)*ht)/((n+1)*ht-n*ht)
+ aux_error = lambda x,t : (exact_rho(x,t)-morley_spacetime(x,t))**2 # L^2 in space
+ grad_morley_spacetime = lambda x,t : grad_morley_interpol1(x)*(t-n*ht)/((n+1)*ht-n*ht) - grad_morley_interpol0(x)*(t-(n+1)*ht)/((n+1)*ht-n*ht)
+ grad_aux_error = lambda x,t : np.linalg.norm(grad_exact_rho(x,t)-grad_morley_spacetime(x,t))**2 # H^1 in space
+
+ # COMPUTE L2 NORM in space
+ aux = lambda x : aux_error(x,n*ht) # the maximum is attaind at time steps t^n = n*ht
+ intL2 = my.get_integral_Omega(aux) # integral over total domain
+ L2_array.append(intL2)
+
+ if n == Nt-1 :
+ aux = lambda x : aux_error(x,(n+1)*ht) # the maximum is attaind at time steps t^n = n*ht
+ intL2 = my.get_integral_Omega(aux) # integral over total domain
+ L2_array.append(intL2)
+
+ # COMPUTE L2H1 NORM
+ maxiter = np.lcm.reduce([1, 2, 3, 4, 6 ,8 , 16, 32, 64]) # get maxiter as lowest common multiple of numbers of time steps used
+ for m in range(int(maxiter/Nt)) : # subintervals
+ I = np.array([(maxiter/Nt*n+m)*T/maxiter,((maxiter/Nt*n+m)+1)*T/maxiter])
+
+ #quadrature rule on each subinterval
+ areaI = I[1]-I[0]
+ trafo = lambda y : y*(I[1]-I[0])/2+(I[1]+I[0])/2 # y in [-1,1]
+ ti = np.array([-np.sqrt(1/3),np.sqrt(1/3)]) #on [-1,1]
+ wi = [1,1]
+ val = 0
+ for k in range(len(wi)) :
+ grad_aux = lambda x : grad_aux_error(x,trafo(ti[k]))
+ intH1 = my.get_integral_Omega(grad_aux)
+ val += wi[k]*intH1
+
+ integral = val/2*areaI # get integral
+
+ L2H1_sum += np.array(integral)
+
+LinfL2 = np.max(L2_array) # compute Linf in time
+
+# SAVE 'EXACT' ERROR IN LinfL2-L2H1 NORM
+data = []
+data.append(L2_array)
+data.append(LinfL2)
+data.append(L2H1_sum)
+print(data)
+pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_exact error.p'
+pickle.dump(data,open('pickle_files/'+pickle_name,'wb')) # store data
+
+progress = 100
+print("%.2f" % progress, 'procent of progress made.')
+elapsed = time.time() - tic
+print('This took',"%.2f" % round(elapsed/60, 2), 'minutes.')
\ No newline at end of file
diff --git a/2d/main_morley.py b/2d/main_morley.py
new file mode 100644
index 0000000000000000000000000000000000000000..929dc8ade91ca481c9afe62491b2161d304ec614
--- /dev/null
+++ b/2d/main_morley.py
@@ -0,0 +1,304 @@
+import myfun as my
+import numpy as np
+import time
+import pickle
+
+# ---------------------------------------------------- FV-FE ALGORITHM -------------------------------------------------------------------------
+# This script performs a FV-FE Algroithm, displayed in the master thesis, for the parabolic-elliptic Keller-Segel system with an initial datum and
+# the homogeneous Neumann boundary condition
+
+tic = time.time()
+
+# SETTINGS FOR SCHEME
+nodal_avrg_choice = 'least-squares' # 'arithmetic-mean' or 'least-squares' -> see nasa paper
+method = 'wasserstein' # using log-mean convective flux
+boundary = 'Neumann'
+interpol = 'morley'
+
+test = 'manufactured1' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
+fineness = 3 # number of refinements of mesh before calculating numerical solution
+Nt = 8 # number of time steps
+
+print('FV-FE Algorithm')
+print(interpol+'_'+boundary+'_'+method+'_'+test)
+print('fineness',fineness)
+print('Nt',Nt)
+
+# initialize different test cases
+if test == 'blow-up' :
+# CONVECTION DEOMINATED REGIME i.e. BLOW-UP
+
+ def rho0(x) :
+ val = 10**(3)*np.exp(-((x[0]-0.5)**2+(x[1]-0.5)**2)/10**(-2))
+ return val
+
+ def f(x,t) :
+ val = 0
+ return val
+
+ def g(x,t) :
+ val = 0
+ return val
+
+ T = 0.0045 # time interval [0,T]
+
+elif test == 'diff' :
+# DIFFUSION DOMINATED REGIME
+ def rho0(x) :
+ val = 1.3*np.exp(-25*(x[0]-1/2)**2-25*(x[1]-1/2)**2)
+ return val
+
+ def f(x,t) :
+ val = 0
+ return val
+
+ def g(x,t) :
+ val = 0
+ return val
+
+ T = 0.05 # time interval [0,T]#
+
+elif test == 'manufactured' :
+# MANUFACTURED SOLUTION WITH T = 0.5
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def delt_exact_rho(x,t) :
+ val = -1.3/(1+t)**2*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def grad_exact_rho(x,t) :
+ val = -1.3*120/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*np.array([x[0]-1/2,x[1]-1/2])
+ return val
+
+ def grad_exact_c(x) :
+ val = -1.3*120*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*np.array([x[0]-1/2,x[1]-1/2])
+ return val
+
+ def Laplacian_exact_rho(x,t) :
+ val = 1/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*(9048-18720*x[0]+18720*x[0]**2-18720*x[1]+18720*x[1]**2)
+ return val
+
+ def Laplacian_exact_c(x) :
+ val = np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*(9048-18720*x[0]+18720*x[0]**2-18720*x[1]+18720*x[1]**2)
+ return val
+
+
+ def f(x,t) :
+ val = delt_exact_rho(x,t) + np.dot(grad_exact_c(x),grad_exact_rho(x,t)) + exact_rho(x,t)*Laplacian_exact_c(x) - Laplacian_exact_rho(x,t)
+ return val
+
+ def g(x,t) :
+ val = exact_c(x) - Laplacian_exact_c(x) - exact_rho(x,t)
+ return val
+
+ T = 0.5 # time interval [0,T]#
+
+elif test == 'manufactured1' :
+# MANUFACTURED SOLUTION WITH T = 3
+
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def delt_exact_rho(x,t) :
+ val = -1.3/(1+t)**2*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def grad_exact_rho(x,t) :
+ val = -1.3*120/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*np.array([x[0]-1/2,x[1]-1/2])
+ return val
+
+ def grad_exact_c(x) :
+ val = -1.3*120*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*np.array([x[0]-1/2,x[1]-1/2])
+ return val
+
+ def Laplacian_exact_rho(x,t) :
+ val = 1/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*(9048-18720*x[0]+18720*x[0]**2-18720*x[1]+18720*x[1]**2)
+ return val
+
+ def Laplacian_exact_c(x) :
+ val = np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)*(9048-18720*x[0]+18720*x[0]**2-18720*x[1]+18720*x[1]**2)
+ return val
+
+
+ def f(x,t) :
+ val = delt_exact_rho(x,t) + np.dot(grad_exact_c(x),grad_exact_rho(x,t)) + exact_rho(x,t)*Laplacian_exact_c(x) - Laplacian_exact_rho(x,t)
+ return val
+
+ def g(x,t) :
+ val = exact_c(x) - Laplacian_exact_c(x) - exact_rho(x,t)
+ return val
+
+ T = 3 # time interval [0,T]#
+
+else :
+ print('Warning: Wrong string input for test.')
+
+toc = time.time()
+
+# SCHEME/PROBLEM PARAMETERS
+hx = 1/2
+ht = (T)/(Nt)
+
+# GET PRIMAL MESH :
+# initialize zerost mesh for FV approximations
+[tri,points,nppoints,numtri] = my.initialmesh()
+K = nppoints[tri.simplices]
+
+# refine zerost mesh subject to fineness
+for i in range(fineness) :
+ hx = 1/2*hx
+ [tri,points,nppoints,numtri] = my.refinemesh(points, K, numtri)
+ K = nppoints[tri.simplices]
+
+#get dual mesh for FE method to approximate chemical density
+[K_dual, tri_dual, points_dual, nppoints_dual, numtri_dual] = my.getdualmesh(points, tri, fineness, K)
+K_dual = np.array(K_dual)
+
+[tri_dual,K_dual] = my.orientation_dualmesh(tri_dual,K_dual) # enforce positive orientation of all triangles in dual mesh
+
+# compute mesh topology for faster find_simplex routine for dual mesh
+[oriented_edges,structure,edge_to_index] = my.get_edges_structure(np.array(K_dual))
+my.init_indexE(numtri_dual)
+
+elapsed = time.time() - toc
+print('Mesh topology computed in ',"%.2f" % round(elapsed/60, 2), 'minutes.')
+toc = time.time()
+
+# print most important information regarding problem/scheme
+print('numtri',numtri)
+print('numtri_dual',numtri_dual)
+print('Nt',Nt)
+print('hx',hx)
+print('ht',ht)
+
+# pre-allocate variables
+rho = np.zeros([Nt+1,numtri])
+cc = np.zeros([Nt+1,len(points)])
+
+toc = time.time()
+
+A = my.assemble_FE_matrix1(tri_dual,K_dual,points_dual)
+elapsed = time.time() - toc
+print('FE matrix assembled in ',"%.2f" % round(elapsed/60, 2), 'minutes.')
+toc = time.time()
+
+# calculate c0 with rho0 as right-hand side
+c0 = my.getc_FE1(rho0,lambda x: g(x,0) ,tri_dual,K_dual,points_dual,A)
+cc[0][:] = c0
+elapsed = time.time() - toc
+print('Initial chemical concentration approximated via FE scheme in ',"%.2f" % round(elapsed/60, 2), 'minutes.')
+toc = time.time()
+
+# get discrete rho0 by evaluating at circumcenters x_K
+for i in range(len(K)) :
+ CC = my.circumcenter(K[i])
+ rho[0][i] = rho0(CC)
+
+# SAVE INITIAL RHO
+data = []
+data.append(rho[0][:])
+data.append(rho0)
+pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(0)+'.p'
+pickle.dump(data,open('pickle_files/'+pickle_name,'wb')) # store data
+
+hht = [0]
+
+#flag = False # for adaptive time stepping
+for n in range(Nt) : # range(Nt) :
+ toc = time.time()
+
+ progress = n/Nt*100
+ print("%.2f" % progress, 'procent of progress made.')
+
+ # SAVE DATA CHEMICAL CONCENTRATION
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_c at time step'+str(n)+'.p'
+ pickle.dump(cc[n][:],open('pickle_files/'+pickle_name,'wb')) # store data
+
+ # GET RHO VIA FINITE VOLUME SCHEME
+ # set v = gradient c
+ v = lambda x: my.grad_cfunc(cc[n][:],tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) # piecewiese constant
+ elapsed = time.time() - toc
+ print('Convection term initialized via v in ',"%.2f" % round(elapsed/60, 2), 'minutes.')
+ toc = time.time()
+
+ # # ADAPTIVE TIME STEPPING
+ # if (T-hht[-1]) < ht : # make sure we do not overshoot time interval [0,T]
+ # ht = np.abs(T-hht[-1])
+ # flag = True
+ # else :
+ # ht = my.get_timestep(v,K,tri)
+
+ # hht.append(hht[-1]+ht) # list of t^n
+ # print('time-step',hht)
+
+ # EXECUTE FINITE VOLUME SCHEME
+ rho[n+1][:] = my.getrho(rho[n][:],lambda x : f(x,n*ht),ht,v,tri,K,numtri)
+ elapsed = time.time() - toc
+ print('Bacterial denisty computed via FV scheme in ', "%.2f" % round(elapsed/60, 2), 'minutes.')
+ toc = time.time()
+
+ # GET COEFFICIENTS FOR MORLEY INTERPOLATION :
+ FKED = my.getinterpolationRHS(tri,K,rho[n+1][:])
+ elapsed = time.time() - toc
+ print('RHS for Morley coefficients calculated in ', "%.2f" % round(elapsed/60, 2), 'minutes.')
+ toc = time.time()
+
+ qq0 = my.getq0(K,nppoints,numtri,tri,rho[n+1][:],nodal_avrg_choice) # linear interpolation
+ elapsed = time.time() - toc
+ print('Linear interpolation q0 computed in ', "%.2f" % round(elapsed/60, 2), 'minutes.')
+ toc = time.time()
+
+ betaKE = my.getbetaE(K,qq0,FKED) # Jacobi_cdotn is first derivative of componentwise linear interpolation of grad c
+ elapsed = time.time() - toc
+ print('Morley coefficients computed in ', "%.2f" % round(elapsed/60, 2), 'minutes.')
+ toc = time.time()
+
+ # GET INTERPOLANT FOR RHO :
+ if test == 'manufactured1' or 'manufactured' :
+ interpol_rho = lambda x : 0
+ else :
+ interpol_rho = lambda x: my.morley_interpol(K,tri,qq0,betaKE,x,indexK=-1) # morley interpolation
+
+ # SAVE BACTERIAL DENSITY RHO
+ data = []
+ data.append(ht) # time step size Delta t^n := t^{n+1} - t^n
+ data.append(rho[n+1][:])
+ data.append(qq0)
+ data.append(betaKE)
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n+1)+'.p'
+ pickle.dump(data,open('pickle_files/'+pickle_name,'wb')) # store data
+
+ # GET CHEMICAL DENISTY
+ cc[n+1][:] = my.getc_FE1(interpol_rho,lambda x : g(x,(n+1)*ht),tri_dual,K_dual,points_dual,A)
+ elapsed = time.time() - toc
+ print('Chemical concentration approximated via FE scheme in ', "%.2f" % round(elapsed/60, 2), 'minutes.')
+
+ if n == Nt-1 :
+ # SAVE DATA CHEMICAL CONCENTRATION
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_c at time step'+str(n+1)+'.p'
+ pickle.dump(cc[n][:],open('pickle_files/'+pickle_name,'wb')) # store data
+
+progress = 100
+print("%.2f" % progress, 'procent of progress made.')
+elapsed = time.time() - tic
+print('This took',"%.2f" % round(elapsed/60, 2), 'minutes.')
+
diff --git a/2d/myfun.py b/2d/myfun.py
new file mode 100644
index 0000000000000000000000000000000000000000..11ae8ac7aa2641c10011c1808f8c3b037eb955f3
--- /dev/null
+++ b/2d/myfun.py
@@ -0,0 +1,1837 @@
+import numpy as np
+from scipy.spatial import Delaunay
+from scipy.sparse import csr_matrix
+from matplotlib import pyplot as plt
+from scipy.sparse.linalg import spsolve
+
+
+
+# GENERAL FUNCTIONS
+
+def initialmesh() :
+# input : None
+# output : a very coarse triangulation (mesh size = 1/2) with property that all inner angles of the triangles are < pi/2
+
+ points = [[0,0],[1,0],[0,1],[1,1]]
+ nppoints = np.array(points)
+
+ points.append(1/2*(nppoints[0]+nppoints[1]))
+ points.append(1/2*(nppoints[0]+nppoints[2]))
+ points.append(1/2*(nppoints[3]+nppoints[1]))
+ points.append(1/2*(nppoints[3]+nppoints[2]))
+ points.append(1/3*(nppoints[0]+nppoints[1]+nppoints[2]))
+ points.append(1/3*(nppoints[3]+nppoints[1]+nppoints[2]))
+ points.append(29/96*nppoints[1]+67/96*nppoints[2]) # where we use 7/24 as the mean of 1/3 und 1/4, to avoid angles of 90°
+ points.append(67/96*nppoints[1]+29/96*nppoints[2])
+
+ nppoints = np.array(points)
+
+ tri = Delaunay(points) # d=2
+ # tri.simplices gives indices that form the triangles
+ # points[tri.simplices] gives vertices of the mesh elements
+ # see documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html
+
+ numtri = 14 # number of triangles in mesh
+
+ return [tri, points, nppoints, numtri]
+
+def refinemesh(points, K, numtri) :
+# input : points that occur in triangulation, triangulation K (array of triangles), numtri number of triangles in triangulation
+# output : delaunay triangulation tri, points, nppoints (=points in format np.array), number of triangles numtri
+ for i in range(len(K)) :
+ newpoints = [1/2*(K[i][0]+K[i][1]), 1/2*(K[i][0]+K[i][2]), 1/2*(K[i][1]+K[i][2])]
+ numtri += 3
+ for j in range(3) :
+ points.append(newpoints[j])
+
+ auxtpl = list(set([tuple(x) for x in points])) # remove duplicates
+ points = [list(ele) for ele in auxtpl]
+ nppoints = np.array(points)
+
+ tri = Delaunay(points)
+ K = nppoints[tri.simplices]
+
+ return [tri, points,nppoints,numtri]
+
+
+def getangles(K):
+# need this for function circumcenter(...)
+
+# input : triangle K
+# output : radian angles of triangle order as the vertices in K, number of triangles in mesh numtri
+
+ A = K[0]
+ B = K[1]
+ C = K[2]
+
+ # Square of lengths be a2, b2, c2
+ a2 = lengthSquare(B, C)
+ b2 = lengthSquare(A, C)
+ c2 = lengthSquare(A, B)
+
+ # length of sides be a, b, c
+ a = np.sqrt(a2);
+ b = np.sqrt(b2);
+ c = np.sqrt(c2);
+
+ # From Cosine law
+ alpha = np.arccos((b2 + c2 - a2)/(2 * b * c));
+ beta = np.arccos((a2 + c2 - b2)/(2 * a * c));
+ gamma = np.arccos((a2 + b2 - c2)/(2 * a * b));
+
+ angles = [alpha, beta, gamma]
+
+ return angles
+
+def lengthSquare(X, Y):
+# need this for function circumcenter(...)
+
+# input : points X and Y
+# output : squared distance between X and Y
+
+ xDiff = X[0] - Y[0]
+ yDiff = X[1] - Y[1]
+
+ return xDiff * xDiff + yDiff * yDiff
+
+def circumcenter(K) :
+# input : triangle K
+# output : circumcenter CC, which can be characterized as the intersection of the perpendicular bisectors of the triangle K
+
+ angles = getangles(K)
+
+ # circumcenter forular from https://testbook.com/maths/circumcenter-of-a-triangle
+ CC = np.array([(K[0][0]*np.sin(2*angles[0])+K[1][0]*np.sin(2*angles[1])+K[2][0]*np.sin(2*angles[2])),
+ (K[0][1]*np.sin(2*angles[0])+K[1][1]*np.sin(2*angles[1])+K[2][1]*np.sin(2*angles[2]))])
+ CC = 1/(np.sin(2*angles[0])+np.sin(2*angles[1])+np.sin(2*angles[2]))*CC
+
+ return CC
+
+
+def are_neighbors(Ki,Li) :
+ # need this for function getdualmesh(...)
+
+ # input: triangles Ki and Li
+ # output: true or false, indicating if Ki and Li are neighbors
+
+ aux = len(list(set([Ki[0],Ki[1],Ki[2],Li[0],Li[1],Li[2]])))
+
+ if aux == 3 :
+ print('Warning: Triangles are the same')
+ flag = False
+
+ elif aux == 4 : # neighbors by edge
+ flag = True
+
+ elif aux == 5 : # neighbors by node
+ flag = False
+
+ else : # no neighbors
+ flag = False
+
+ return flag
+
+def linear_search(list, x):
+ # need this for function getdualmesh(...)
+ # basic line search algorithm for 2 dim vectors
+ for i in range(len(list)):
+ if list[i][0] == x[0] and list[i][1] == x[1]:
+ return i
+ return -1
+
+def linear_search_edges(list, x):
+ # need this for function find_simlex1(...)
+ # basic line search algorithm for edges
+ for i in range(len(list)):
+ if list[i][0][0] == x[0][0] and list[i][0][1] == x[0][1] and list[i][1][0] == x[1][0] and list[i][1][1] == x[1][1]:
+ return i
+ return -1
+
+def getdualmesh(points, tri, fineness, K):
+ # input: primal mesh and some of its parameters
+ # output: dual mesh for the FE scheme to induce no jumps for (primal) normal derivatives
+
+ points_primal = len(points)
+
+ hx = 1/2**(fineness+1)
+ points_dual = points
+ nppoints = np.array(points) # important that we initialite this here as list points is changed below...
+
+ for i in range(len(K)) :
+ points_dual.append(circumcenter(K[i]))
+
+ nppoints_dual = np.array(points_dual)
+
+ pt_i = -1
+ K_dual = []
+
+ class myclass :
+ def __init__(self) :
+ self.simplices = []
+
+ tri_dual = myclass()
+
+ for pt in nppoints : # go through all points in primal mesh
+ pt_i +=1
+
+ neighbors = find_neighbors(pt_i,tri)
+ taken_care_of = []
+
+ for i in neighbors :
+ taken_care_of.append(i)
+
+ for j in neighbors :
+
+ if j in taken_care_of :
+ continue
+
+ if are_neighbors(tri.simplices[i],tri.simplices[j]) :
+ CCi = circumcenter(K[i])
+ CCj = circumcenter(K[j])
+
+ K_dual.append([pt,CCi,CCj])
+ tri_dual.simplices.append([pt_i,points_primal+i,points_primal+j])
+
+ # extra triangles on boundary
+ if pt[0] == 0 : # bdr edge x = 0
+ if pt[1] == 0: # corner [0,0]
+ x = pt+np.array([hx,0])
+ y = pt+np.array([0,hx])
+
+ midx = 1/2*(pt+x)
+ midy = 1/2*(pt+y)
+
+ i = tri.find_simplex(midx+10**-8*np.array([0,1]))
+ j = tri.find_simplex(midy+10**-8*np.array([1,0]))
+
+ CCi = circumcenter(K[i])
+ CCj = circumcenter(K[j])
+
+ K_dual.append([pt,x,CCi])
+ x_i = linear_search(nppoints,x)
+ tri_dual.simplices.append([pt_i,x_i,points_primal+i])
+
+ K_dual.append([pt,y,CCj])
+ y_i = linear_search(nppoints,y)
+ tri_dual.simplices.append([pt_i,y_i,points_primal+j])
+
+ elif pt[1] == 1 : # corner [0,1]
+ x = pt+np.array([hx,0])
+
+ midx = 1/2*(pt+x)
+
+ i = tri.find_simplex(midx-10**-8*np.array([0,1]))
+
+ CCi = circumcenter(K[i])
+ x_i = linear_search(nppoints,x)
+ tri_dual.simplices.append([pt_i,x_i,points_primal+i])
+ K_dual.append([pt,x,CCi])
+
+ else : # no corner
+ x = pt+np.array([0,hx])
+
+ midx = 1/2*(pt+x)
+
+ i = tri.find_simplex(midx+10**-8*np.array([1,0]))
+
+ CCi = circumcenter(K[i])
+ x_i = linear_search(nppoints,x)
+ tri_dual.simplices.append([pt_i,x_i,points_primal+i])
+ K_dual.append([pt,x,CCi])
+
+ elif pt[0] == 1 : # bdr edge x = 1
+
+ if pt[1] == 0: # corner [1,0]
+ y = pt+np.array([0,hx])
+
+ midy = 1/2*(pt+y)
+
+ j = tri.find_simplex(midy-10**-8*np.array([1,0]))
+
+ CCj = circumcenter(K[j])
+ y_i = linear_search(nppoints,y)
+ tri_dual.simplices.append([pt_i,y_i,points_primal+j])
+ K_dual.append([pt,y,CCj])
+
+ elif pt[1] == 1 : # corner [1,1]
+ nothing = 1
+
+ else : # no corner
+ x = pt+np.array([0,hx])
+
+ midx = 1/2*(pt+x)
+
+ i = tri.find_simplex(midx-10**-8*np.array([1,0]))
+
+ CCi = circumcenter(K[i])
+ x_i = linear_search(nppoints,x)
+ tri_dual.simplices.append([pt_i,x_i,points_primal+i])
+ K_dual.append([pt,x,CCi])
+
+ elif pt[1] == 0: # bdr edge y = 0
+ x = pt+np.array([hx,0])
+
+ midx = 1/2*(pt+x)
+
+ i = tri.find_simplex(midx+10**-8*np.array([0,1]))
+
+ CCi = circumcenter(K[i])
+ x_i = linear_search(nppoints,x)
+ tri_dual.simplices.append([pt_i,x_i,points_primal+i])
+ K_dual.append([pt,x,CCi])
+
+ elif pt[1] == 1: # bdr edge y = 1
+ x = pt+np.array([hx,0])
+
+ midx = 1/2*(pt+x)
+
+ i = tri.find_simplex(midx-10**-8*np.array([0,1]))
+
+ CCi = circumcenter(K[i])
+ x_i = linear_search(nppoints,x)
+ tri_dual.simplices.append([pt_i,x_i,points_primal+i])
+ K_dual.append([pt,x,CCi])
+
+
+ numtri_dual = len(K_dual)
+
+ return [K_dual, tri_dual, points_dual, nppoints_dual, numtri_dual]
+
+def orientation_dualmesh(tri,K) :
+# input: dual mesh
+# output: dual mesh, s.t. all triangles are oriented positively, i.e. the vertices that define a triangle are ordered counterclockwise
+
+ poss = [[0,1,2],[0,2,1],[1,0,2],[1,2,0],[2,1,0],[2,0,1]]
+ for i in range(len(K)) :
+ for k in range(12) :
+ j = np.mod(k,6)
+ flag = -2
+
+ tau1 = K[i][poss[j][1]]-K[i][poss[j][0]]
+ tau2 = K[i][poss[j][2]]-K[i][poss[j][0]]
+
+ if (tau1[0] > 0 and tau1[1] > 0) :
+ flag += 1
+
+ if k < 6 :
+ tau1_perp1 = [-tau1[1],tau1[0]]
+ else :
+ tau1_perp1 = [tau1[1],-tau1[0]]
+
+ if np.dot(tau1_perp1,tau2) > 0 :
+ flag += 1
+
+ if flag == 0 :
+
+ if k < 6 :
+ oldtri = np.array(tri.simplices[i])
+ K[i] = K[i][poss[j]]
+ tri.simplices[i] = oldtri[poss[j]]
+ else :
+ aux_poss = [poss[j][0],poss[j][2],poss[j][1]]
+ oldtri = np.array(tri.simplices[i])
+ K[i] = K[i][aux_poss]
+ tri.simplices[i] = oldtri[aux_poss]
+
+ break
+
+ elif k == 11 :
+ print('Warning: orientation does not work')
+
+ return [tri,K]
+
+
+def RightOf(x,E) :
+ # need tis for function get_edges_structure(...)
+
+ # says if a point is on the right of an oriented edge
+ val = (E[1][0] - E[0][0])*(x[1] - E[0][1]) - (E[1][1] - E[0][1])*(x[0] - E[0][0])
+ if val < 0 :
+ return True
+ else :
+ return False
+
+def get_edges_structure(K) :
+# need this for more efficient find_simplex1 for the dual mesh
+
+# input: dual mesh in terms of K
+# output: mesh topology of the dual mesh
+
+ oriented_edges = []
+ structure = []
+ edge_to_index = []
+
+ for i in range(len(K)) :
+ indices = [[1,2],[0,2],[0,1]]
+ for j in range(3) :
+
+ E = K[i][indices[j]]
+ pt = K[i][j]
+
+ if RightOf(pt,E) :
+ E = [E[1],E[0]]
+
+ oriented_edges.append(E)
+ structure.append([[E[0],pt],[pt,E[1]]])
+ edge_to_index.append(i)
+
+ return [oriented_edges,structure,edge_to_index]
+
+
+def find_neighbors(pindex, tri):
+# input : index of point pindex, delaunay triangulation tri
+# output : indices of triangles that point pindex is part of
+
+ #neighborstri = []
+ neighborstri_index = []
+ i = -1
+ for simplex in tri.simplices:
+ i += 1
+ if pindex in simplex:
+ #neighborstri.append(simplex)
+ neighborstri_index.append(i)
+
+ return neighborstri_index
+
+
+def min_dist(E,x) :
+# need this for function find_simplex1(...)
+ # get minimal distance between edge E and point x
+ val = np.abs(np.cross(E[1]-E[0], E[0]-x))/np.linalg.norm(E[1]-E[0])
+ return val
+
+def find_simplex1(x,oriented_edges,structure,indexE) :
+#input: point x, mesh topology of the dual mesh
+# output: edge E on which x lies or x lies in the triangle left of E
+
+ # algorithm cannot deal with edges so we weak it a bit
+ if x[0] == 0 :
+ x[0] = 10**-10
+ if x[1] == 0 :
+ x[1] = 10**-10
+ if x[0] == 1 :
+ x[0] = 1-10**-10
+ if x[1] == 1 :
+ x[1] = 1-10**-10
+
+ E = oriented_edges[indexE]
+
+ if RightOf(x, E) :
+ E = [E[1],E[0]]
+
+ indexE = linear_search_edges(oriented_edges,E)
+
+ list_indices = []
+
+ while True :
+ if indexE in list_indices :
+ print('x',x)
+ print('Warning: loop.')
+ break
+ list_indices.append(indexE)
+
+ if (x[0] == E[0][0] and x[1] == E[0][1]) or (x[0] == E[1][0] and x[1] == E[1][1]) :
+ return indexE
+ else :
+ whichop = 0
+ if not(RightOf(x, structure[indexE][0])) : # Onext
+ whichop += 1
+ if not(RightOf(x, structure[indexE][1])) : #Dprev
+ whichop += 2
+
+ if whichop == 0 :
+ return indexE
+ elif whichop == 1:
+ E = structure[indexE][0] # Onext
+ elif whichop == 2:
+ E = structure[indexE][1] # Dprev
+ elif whichop == 3:
+ if min_dist(structure[indexE][0], x) < min_dist(structure[indexE][1], x) :
+ E = structure[indexE][0] # Onext
+ else :
+ E = structure[indexE][1] # Dprev
+
+ indexE = linear_search_edges(oriented_edges,E)
+
+def init_indexE(numtri_dual) :
+# initialize global variable indexE for more performent find_simplex1
+ global indexE
+ indexE = round(3*numtri_dual/2)
+ # global counter
+ # counter = np.zeros(3*numtri_dual)
+ print('indexE',indexE)
+
+
+def unitouternormal(K,E) :
+# input: element K and edge E in terms of vertices, 2-dim vector x that lies on E
+# output: unit outward normal vector to K along E
+
+ # calculate normal
+ edge = E[1]-E[0]
+ normal = np.array([edge[1],-edge[0]]) # np.dot(normal,edge) =!= 0
+ normal = 1/np.linalg.norm(normal)*normal # normalization
+
+ # check for orientation
+ center = (K[0]+K[1]+K[2])/3 # convex compination of vertices (always in interior triangle)
+ lengthplus = np.linalg.norm((E[0]+E[1])/2+normal-center)
+ lengthminus = np.linalg.norm((E[0]+E[1])/2-normal-center)
+
+ if lengthplus<lengthminus :
+ normal = -normal
+
+ return normal
+
+
+def tritrafo(K,L,pt) :
+# need this for integrating via function avrgK(...)
+
+# input : triangles K and L, point pt from K
+# output : pt transformed to L
+
+ A = np.array([[K[0][0], K[1][0], K[2][0]], # | xa1 xa2 xa3 |
+ [K[0][1], K[1][1], K[2][1]], #A =| ya1 ya2 ya3 |
+ [ 1, 1, 1]]) # | 1 1 1 |
+ B = np.array([[L[0][0], L[1][0], L[2][0]], # | xa1 xa2 xa3 |
+ [L[0][1], L[1][1], L[2][1]], #A =| ya1 ya2 ya3 |
+ [ 1, 1, 1]]) # | 1 1 1 |
+
+
+ invA = np.linalg.solve(A,np.eye(3))
+ #invA = np.linalg.inv(A)
+ M = np.matmul(B,invA)
+
+ auxx = np.array([pt[0],pt[1],1])
+
+ auxtrafo = np.matmul(M,auxx)
+ trafo = [auxtrafo[0],auxtrafo[1]]
+
+ return [trafo,M]
+
+def avrgK(K,g) :
+# input : element K in terms of vertices, function g
+# output : average of g:IR²->IR over K, area of K
+ S = np.cross(K[1]-K[0],K[2]-K[0]) # B = K[0], C = K[1], A = K[2]
+ areaK = np.abs(S)/2 # get area of triangle K
+
+ # wi = [1] (exactness 1 -> Order 2)
+ # xi = np.array([[1/3,1/3]])
+
+ # from https://mathsfromnothing.au/triangle-quadrature-rules/?i=1
+ # wi = [0.223381589678011,0.223381589678011,0.223381589678011,0.109951743655322,0.109951743655322,0.109951743655322]
+ # xi = np.array([[0.445948490915965,0.108103018168070],
+ # [0.445948490915965,0.445948490915965],
+ # [0.108103018168070,0.445948490915965],
+ # [0.091576213509771,0.816847572980459],
+ # [0.091576213509771,0.091576213509771],
+ # [0.816847572980459,0.091576213509771]])
+
+
+ wi = [1/6,1/6,1/6] # weights*area_unitK(=1/2) (Exactness 2 -> Order 3)
+ xi = np.array([[1/6,2/3],[1/6,1/6],[2/3,1/6]]) # points in unit triangle
+
+ unitK = [[0,0],[1,0],[0,1]]
+
+ val = 0
+ for i in range(len(wi)) :
+ [trafoxi,M] = tritrafo(unitK,K,xi[i])
+ val += wi[i]*g(trafoxi)
+
+ avrg = val*2
+ integral = avrg*areaK
+
+ return [integral,avrg,areaK]
+
+def avrgE(E,g) :
+ # input : edge E in terms of vertices, function g
+ # output : average of g:IR²->IR over E, area of E
+
+ areaE = np.linalg.norm(E[0]-E[1])
+
+ trafo = lambda y : y*(E[1]-E[0])/2+(E[1]+E[0])/2 # y in [-1,1]
+
+ xi = np.array([-np.sqrt(1/3),np.sqrt(1/3)]) #on [-1,1]
+ wi = [1,1]
+
+ # xi = np.array([-np.sqrt(3/5),0,np.sqrt(3/5)])
+ # wi = [5/9,8/9,5/9]
+
+ val = 0
+ for k in range(len(wi)) :
+ val += wi[k]*g(trafo(xi[k])) # get integral
+
+ avrg = val/2 # get average
+ integral = avrg*areaE
+
+ return [integral,avrg,areaE]
+
+def quadrature(I,g) :
+# inupt: interval I = [a,b], function g
+# output: approximated integral of g over I
+
+ areaI = I[1]-I[0]
+
+ trafo = lambda y : y*(I[1]-I[0])/2+(I[1]+I[0])/2 # y in [-1,1]
+
+ xi = np.array([-np.sqrt(1/3),np.sqrt(1/3)]) #on [-1,1]
+ wi = [1,1]
+
+ # xi = [0]
+ # wi = [2]
+
+ # xi = np.array([-np.sqrt(3/5),0,np.sqrt(3/5)]) #on [-1,1]
+ # wi = [5/9,8/9,5/9]
+
+ val = 0
+ for k in range(len(wi)) :
+ val += wi[k]*g(trafo(xi[k])) # get integral
+
+ avrg = val/2 # get average
+ integral = avrg*areaI
+
+ return integral
+
+
+def get_integral_Omega(g) :
+# input: function g
+# output: integral of g over Omega
+
+ # get initial mesh for quadrature (different to one that we use to calculate numerical solution to avoid unwanted effects due to projection errors)
+ points = [[0,0],[1,0],[0.5,0.5],[0,1],[1,1]]
+ nppoints = np.array(points)
+ tri = Delaunay(points)
+ K = nppoints[tri.simplices]
+ numtri = 2
+
+ # refine the mesh until it is finer than mesh used to calculate numerical solution, here 5
+ for i in range(5) :
+ [tri,points,nppoints,numtri] = refinemesh(points, K, numtri)
+ K = nppoints[tri.simplices]
+
+ # calculate integral over Omega as sum over inetrgals on K
+ integral = 0
+ for i in range(numtri) :
+ [aux,avrg,areaK] = avrgK(K[i],g)
+ integral += aux
+
+ return integral
+
+
+def newton_method(f,grad_f,x0,tol,maxiter) :
+# classical newton method
+# Input: function f, Jacobian of f in sparse form, starting point x0, tolerance tol and maximal amount of iteration maxiter
+# output: root of f
+
+ xnew = x0
+ for i in range(maxiter) :
+ xold = xnew
+
+ delx = spsolve(grad_f(xold), -f(xold), permc_spec=None, use_umfpack=True)
+
+ xnew = xold + delx
+
+ if np.linalg.norm(xold-xnew) < tol :
+ break
+
+ print('Newton steps:',i)
+
+ if i == maxiter-1 :
+ print('Warning: Newton-method reached maxiter.')
+
+ return xnew
+
+
+def bubbleK(K,pt) :
+# input : triangle K, point pt
+# output : element bubble function
+
+ # get reference bubble functions on unit triangle
+ refb1 = lambda x: 1-x[0]-x[1]
+ refb2 = lambda x: x[0]
+ refb3 = lambda x: x[1]
+ refb = lambda x: refb1(x)*refb2(x)*refb3(x)*27
+
+ # transform the triangle to get general bubble function
+ unitK = [[0,0],[1,0],[0,1]]
+ [trafo,M] = tritrafo(K,unitK,pt)
+
+ val = refb(trafo)
+
+ return val
+
+def gradbubbleK(K,pt) :
+# input : triangle K, point x
+# output : gradient of element bubble function
+
+ # get reference bubble functions on unit triangle
+ refb1 = lambda x: 1-x[0]-x[1]
+ refb2 = lambda x: x[0]
+ refb3 = lambda x: x[1]
+ gradrefb1 = lambda x: [-1,-1]
+ gradrefb2 = lambda x: [1,0]
+ gradrefb3 = lambda x: [0,1]
+ gradrefb = lambda x: (np.multiply(gradrefb1(x),refb2(x))*refb3(x)+refb1(x)*np.multiply(gradrefb2(x),refb3(x))+refb1(x)*np.multiply(gradrefb3(x),refb2(x)))*27
+
+ # transform the triangle to get general bubble function
+ unitK = [[0,0],[1,0],[0,1]]
+ [trafo,M] = tritrafo(K,unitK,pt)
+
+ # chain rule
+ DM = np.array([[M[0,0],M[0,1]],[M[1,0],M[1,1]]]) # = M[0:1,0:1]
+ val = np.matmul(gradrefb(trafo),DM) # order of multiplikation important!
+
+ return val
+
+def LaplacebubbleK(K,pt) :
+# input : triangle K, point x
+# output : Laplace of element bubble function
+
+ # get reference bubble functions on unit triangle
+ refb1 = lambda x: 1-x[0]-x[1]
+ refb2 = lambda x: x[0]
+ refb3 = lambda x: x[1]
+ gradrefb1 = lambda x: [-1,-1]
+ gradrefb2 = lambda x: [1,0]
+ gradrefb3 = lambda x: [0,1]
+ Hesse_refb = lambda x: [[27*(refb2(x)*gradrefb1(x)[0]*gradrefb3(x)[0]+refb3(x)*gradrefb1(x)[0]*gradrefb2(x)[0]+refb1(x)*gradrefb3(x)[0]*gradrefb2(x)[0]+refb2(x)*gradrefb3(x)[0]*gradrefb1(x)[0]+refb1(x)*gradrefb2(x)[0]*gradrefb3(x)[0]+refb3(x)*gradrefb2(x)[0]*gradrefb1(x)[0]) , 27*(refb2(x)*gradrefb1(x)[1]*gradrefb3(x)[0]+refb3(x)*gradrefb1(x)[1]*gradrefb2(x)[0]+refb1(x)*gradrefb3(x)[1]*gradrefb2(x)[0]+refb2(x)*gradrefb3(x)[1]*gradrefb1(x)[0]+refb1(x)*gradrefb2(x)[1]*gradrefb3(x)[0]+refb3(x)*gradrefb2(x)[1]*gradrefb1(x)[0])],
+ [27*(refb2(x)*gradrefb1(x)[0]*gradrefb3(x)[1]+refb3(x)*gradrefb1(x)[0]*gradrefb2(x)[1]+refb1(x)*gradrefb3(x)[0]*gradrefb2(x)[1]+refb2(x)*gradrefb3(x)[0]*gradrefb1(x)[1]+refb1(x)*gradrefb2(x)[0]*gradrefb3(x)[1]+refb3(x)*gradrefb2(x)[0]*gradrefb1(x)[1]) , 27*(refb2(x)*gradrefb1(x)[1]*gradrefb3(x)[1]+refb3(x)*gradrefb1(x)[1]*gradrefb2(x)[1]+refb1(x)*gradrefb3(x)[1]*gradrefb2(x)[1]+refb2(x)*gradrefb3(x)[1]*gradrefb1(x)[1]+refb1(x)*gradrefb2(x)[1]*gradrefb3(x)[1]+refb3(x)*gradrefb2(x)[1]*gradrefb1(x)[1])]]
+
+ # transform the triangle to get general bubble function
+ unitK = [[0,0],[1,0],[0,1]]
+ [trafo,M] = tritrafo(K,unitK,pt)
+
+ # double chain rule
+ DM1 = np.array([M[0,0],M[1,0]])
+ DM2 = np.array([M[0,1],M[1,1]])
+ val = np.dot(DM1,np.matmul(Hesse_refb(trafo),DM1)) + np.dot(DM2,np.matmul(Hesse_refb(trafo),DM2)) # due to Laplace chain rule
+
+ return val
+
+
+def bubbleE(K,index_notE,pt) :
+# input : triangle K, indexKpt gives index of point of K that is not part of E (therefore uniquely defining E), point x
+# output :edge bubble function bE
+
+ # get reference bubble functions on unit triangle
+ refb1 = lambda x: 1-x[0]-x[1]
+ refb2 = lambda x: x[0]
+ refb3 = lambda x: x[1]
+
+ if index_notE == 0 : # remove function that is associated to point not belonging to the edge
+ refb = lambda x: refb2(x)*refb3(x)*4 # eher 4
+ elif index_notE == 1 :
+ refb = lambda x : refb1(x)*refb3(x)*4
+ else :
+ refb = lambda x: refb1(x)*refb2(x)*4
+
+ # transform the triangle to get general bubble function
+ unitK = [[0,0],[1,0],[0,1]]
+ [trafo,M] = tritrafo(K,unitK,pt)
+
+ return refb(trafo)
+
+def gradbubbleE(K,index_notE,pt) :
+ # input : triangle K, indexKpt gives index of point of K that is not part of E (therefore uniquely defining E), point x
+ # output : gradient of edge bubble function bE
+
+ # get reference bubble functions on unit triangle
+ refb1 = lambda x: 1-x[0]-x[1]
+ refb2 = lambda x: x[0]
+ refb3 = lambda x: x[1]
+ gradrefb1 = lambda x: [-1,-1]
+ gradrefb2 = lambda x: [1,0]
+ gradrefb3 = lambda x: [0,1]
+
+ if index_notE == 0 : # remove function that is associated to point not belonging to the edge
+ gradrefb = lambda x: (np.multiply(refb2(x),gradrefb3(x))+np.multiply(gradrefb2(x),refb3(x)))*4
+ elif index_notE == 1 :
+ gradrefb = lambda x : (np.multiply(refb1(x),gradrefb3(x))+np.multiply(gradrefb1(x),refb3(x)))*4
+ else :
+ gradrefb = lambda x: (np.multiply(refb1(x),gradrefb2(x))+np.multiply(gradrefb1(x),refb2(x)))*4
+
+ # transform the triangle to get general bubble function
+ unitK = [[0,0],[1,0],[0,1]]
+ [trafo,M] = tritrafo(K,unitK,pt)
+
+ # chain rule
+ DM = [[M[0,0],M[0,1]],[M[1,0],M[1,1]]] # = M[0:1,0:1]
+ val = np.matmul(gradrefb(trafo),DM) # order of multiplikation important!
+
+ return val
+
+def LaplacebubbleE(K,index_notE,pt) :
+ # input : triangle K, indexKpt gives index of point of K that is not part of E (therefore uniquely defining E), point x
+ # output : Laplace of edge bubble function bE
+
+ # get reference bubble functions on unit triangle
+ gradrefb1 = lambda x: [-1,-1]
+ gradrefb2 = lambda x: [1,0]
+ gradrefb3 = lambda x: [0,1]
+
+ if index_notE == 0 : # remove function that is associated to point not belonging to the edge
+ Hesse_refb = lambda x: [[4*(2*gradrefb2(x)[0]*gradrefb3(x)[0]) , 4*(gradrefb2(x)[0]*gradrefb3(x)[1]+gradrefb2(x)[1]*gradrefb3(x)[0])],
+ [4*(gradrefb2(x)[1]*gradrefb3(x)[0]+gradrefb2(x)[0]*gradrefb3(x)[1]) , 4*(2*gradrefb2(x)[1]*gradrefb3(x)[1])]]
+ elif index_notE == 1 :
+ Hesse_refb = lambda x: [[4*(2*gradrefb1(x)[0]*gradrefb3(x)[0]) , 4*(gradrefb1(x)[0]*gradrefb3(x)[1]+gradrefb1(x)[1]*gradrefb3(x)[0])],
+ [4*(gradrefb1(x)[1]*gradrefb3(x)[0]+gradrefb1(x)[0]*gradrefb3(x)[1]) , 4*(2*gradrefb1(x)[1]*gradrefb3(x)[1])]]
+ else :
+ Hesse_refb = lambda x: [[4*(2*gradrefb1(x)[0]*gradrefb2(x)[0]) , 4*(gradrefb1(x)[0]*gradrefb2(x)[1]+gradrefb1(x)[1]*gradrefb2(x)[0])],
+ [4*(gradrefb1(x)[1]*gradrefb2(x)[0]+gradrefb1(x)[0]*gradrefb2(x)[1]) , 4*(2*gradrefb1(x)[1]*gradrefb2(x)[1])]]
+
+ # transform the triangle to get general bubble function
+ unitK = [[0,0],[1,0],[0,1]]
+ [trafo,M] = tritrafo(K,unitK,pt)
+
+ # double chain rule
+ DM1 = np.array([M[0,0],M[1,0]])
+ DM2 = np.array([M[0,1],M[1,1]])
+ val = np.dot(DM1,np.matmul(Hesse_refb(trafo),DM1)) + np.dot(DM2,np.matmul(Hesse_refb(trafo),DM2)) # due to Laplace chain rule
+ return val
+
+
+def get_timestep(v,K,tri) :
+ # input : v = grad_c, triangulation
+ # output: time step
+
+ # motivated by positivity proof / CFL condition in Kolbe
+ # go through all CCs, find adjacent triangles in dual_mesh, calculate grad_c*n_E on each of these adjacent triangles
+ # not used in simulations
+
+ upper_list = []
+
+ for i in range(len(K)) :
+
+ neighbors = tri.neighbors[i] # get neighbors of triangle K[i]
+
+ indices = [[1,2],[0,2],[0,1]]
+ val = 0
+
+ [integral,avrg,areaK] = avrgK(K[i],lambda x: np.dot(x,x))
+
+ for j in range(3):
+
+ if neighbors[j] >= 0 : # i.e. edge E is interior
+
+ E = K[i][indices[j]]
+ midE = 1/2*(E[0]+E[1])
+ normalKE = unitouternormal(K[i],E)
+
+ [integral,avrg,areaE] = avrgE(E,lambda x: np.dot(x,x))
+
+ val =+ np.abs(np.dot(v(midE),normalKE))*areaE
+
+ else : # i.e. edge E is part of the boundary -> vKE = 0 for Neumann bdr condition
+
+ val =+ 0
+
+ upper_list.append(areaK/val)
+
+
+ upper = np.min(upper_list)
+
+ return upper
+
+
+
+# CHEMICAL DENSITY RELATED FUNCTIONS
+
+def hat(i,tri,points,x,oriented_edges,structure,edge_to_index,indexK) :
+# 2D hat function on triangle
+#input: index i of vertex in points, triangulation tri,K , point x, indexK: wenn -1 muss noch gesucht werden wenn > -1 dann index von triangle
+#output: value of hat function at x
+
+ # get triangle in which is x
+ if indexK < -0.5 :
+ global indexE
+ indexE = find_simplex1(x,oriented_edges,structure,indexE)
+ k = edge_to_index[indexE]
+ else :
+ k = indexK
+
+ if i in tri.simplices[k] :
+
+ # get adjacent points of i
+ L = [points[j] for j in tri.simplices[k] if not j==i ]
+
+ M = np.zeros([3,3])
+ rhs = np.zeros(3)
+ rhs[0] = 1 # nodal values for hat functions
+
+ # get system of equations
+ M[0][0] = points[i][0] # a*x+b*y +c = y
+ M[0][1] = points[i][1]
+ M[0][2] = 1
+
+ M[1][0] = L[0][0]
+ M[1][1] = L[0][1]
+ M[1][2] = 1
+
+ M[2][0] = L[1][0]
+ M[2][1] = L[1][1]
+ M[2][2] = 1
+
+ coeff = np.linalg.solve(M,rhs) # solve linear system of equations
+
+ y = coeff[0]*x[0]+coeff[1]*x[1]+coeff[2]
+
+ else :
+ y = 0 # zero if outside of stencil
+
+ return y
+
+def dxhat(i,tri,points,x,oriented_edges,structure,edge_to_index,indexK) :
+# gradient 2d hat function on triangle
+#input: index i of vertex in points, triangulation tri,K , point x, indexK, wenn -1 muss noch gesucht werden wenn > -1 dann index von triangle
+#output: value of gradient of hat function at x
+
+ if indexK < -0.5 :
+ global indexE
+ indexE = find_simplex1(x,oriented_edges,structure,indexE)
+ k = edge_to_index[indexE]
+ else :
+ k = indexK
+
+ if i in tri.simplices[k] :
+
+ # get adjacent points of i
+ L = [points[j] for j in tri.simplices[k] if not j==i]
+
+
+ M = np.zeros([3,3])
+ rhs = np.zeros(3)
+ rhs[0] = 1 # nodal values for hat functions
+
+ # get system of equations
+ M[0][0] = points[i][0] # a*x+b*y +c = y
+ M[0][1] = points[i][1]
+ M[0][2] = 1
+
+ M[1][0] = L[0][0]
+ M[1][1] = L[0][1]
+ M[1][2] = 1
+
+ M[2][0] = L[1][0]
+ M[2][1] = L[1][1]
+ M[2][2] = 1
+
+ coeff = np.linalg.solve(M,rhs) # solve linear system of equations
+
+ grad = np.array([coeff[0],coeff[1]]) # gradient
+
+ else :
+ grad = np.zeros(2) # zero if x is outside of stencil
+
+ return grad
+
+
+def assemble_FE_matrix1(tri_dual,K_dual,points_dual):
+# see Fig. 3.14. in Bartels: Numerical Approximation of PDEs, 2015
+# input : triangulation
+# output : (affine linear) FE matrix
+
+ # initialize objects
+ iter = 0
+ iter_max = 9*len(K_dual)
+
+ I = np.zeros(iter_max)
+ J = np.zeros(iter_max)
+ X_diff = np.zeros(iter_max)
+ X_reac = np.zeros(iter_max)
+
+ m_loc = 1/12*(np.ones([3,3]) + np.array([[1,0,0],[0,1,0],[0,0,1]])) # formular from Bartels book
+
+ for i in range(len(K_dual)) :
+
+ X_K = np.array([[1, 1, 1],[K_dual[i][0][0], K_dual[i][1][0], K_dual[i][2][0]],[K_dual[i][0][1], K_dual[i][1][1], K_dual[i][2][1]]])
+ rhs = np.array([[0,0],[1,0],[0,1]])
+ grads_K = np.linalg.solve(X_K,rhs) # compute gradients
+ vol_K = np.linalg.det(X_K)/2 # compute areaK
+
+ # get matrix
+ for m in [0,1,2] :
+ for n in [0,1,2] :
+ I[iter] = tri_dual.simplices[i][m]
+ J[iter] = tri_dual.simplices[i][n]
+ X_diff[iter] = vol_K * np.dot(grads_K[m],grads_K[n]) # diffusion contributions (stiffness matrix)
+ X_reac[iter] = vol_K*m_loc[m][n] # reaction conrtibutions (mass matrix)
+ iter += 1
+
+ sparseA_diff = csr_matrix((X_diff[:iter], (I[:iter], J[:iter])), shape=(len(points_dual),len(points_dual)))
+
+ sparseA_reac = csr_matrix((X_reac[:iter], (I[:iter], J[:iter])), shape=(len(points_dual),len(points_dual)))
+
+ sparseA = sparseA_diff + sparseA_reac
+
+ return sparseA
+
+def getc_FE1(interpol_rho,g,tri_dual,K_dual,points_dual,sparseA) :
+# input : interpolation of FV solution rho, rhs g for manufactured solutions, dual triangulation (tri,K,points) and the FE matrix
+# output: coefficients FE solution c
+
+ b = np.zeros(len(points_dual))
+
+ aux = lambda x : interpol_rho(x)+g(x)
+
+ for i in range(len(K_dual)) :
+
+ # get objects from Bartels book
+ X_K = np.array([[1, 1, 1],[K_dual[i][0][0], K_dual[i][1][0], K_dual[i][2][0]],[K_dual[i][0][1], K_dual[i][1][1], K_dual[i][2][1]]])
+ vol_K = np.linalg.det(X_K)/2
+ mid_K = 1/3*(K_dual[i][0]+K_dual[i][1]+K_dual[i][2]) # midpoint of triangle K_dual
+
+ for m in [0,1,2] :
+ b[tri_dual.simplices[i][m]] += 1/3* vol_K * aux(mid_K) # using midpoint rule
+
+ c = spsolve(sparseA, b, permc_spec=None, use_umfpack=True) # solve system of equations
+
+ # Then FE solution is c_h^n(x) = sum_i (c_i*hat_i(x))
+ return c
+
+
+def cfunc(c,tri,K,points,x,oriented_edges,structure,edge_to_index) :
+# input: coefficients FE sol., point x, dual triangulation (tri,K,points), objects to make find_simplex for dual meshes more efficient
+# output: function value of FE solution: c_h^n(x) = sum_i (c_i*hat_i(x))
+
+ val = 0
+ global indexE
+ indexE = find_simplex1(x,oriented_edges,structure,indexE)
+ k = edge_to_index[indexE] # get index of triangle in which x lives
+
+ for i in range(len(c)) : # sum over all coefficients c (stencil would be enough)
+ val += c[i]*hat(i,tri,points,x,oriented_edges,structure,edge_to_index,indexK = k)
+
+ return val
+
+def grad_cfunc(c,tri,K,points,x,oriented_edges,structure,edge_to_index) :
+# input: coefficients FE sol., point x, dual triangulation (tri,K,points), objects to make find_simplex for dual meshes more efficient
+# output: gradient of FE solution: c_h^n(x) = sum_i (c_i*hat_i(x))
+
+ val = 0
+ global indexE
+ indexE = find_simplex1(x,oriented_edges,structure,indexE)
+ k = edge_to_index[indexE]
+
+ for i in range(len(c)) : # sum over all coefficients c (stencil would be enough)
+ val += c[i]*dxhat(i,tri,points,x,oriented_edges,structure,edge_to_index,indexK = k)
+ return val
+
+
+def get_c_plot(pltx,tri,K,points,n,ht,cc,clim,fineness,interpol,oriented_edges,structure,edge_to_index) :
+# input: plot resolution pltx, triangulation, time step n and step size ht, FE sol. cc, fineness, interpolation type and objects to make find_simplex for the dual mesh more effient
+# output: plot of FE solution function c
+
+ # discretize space into pixels
+ x = np.linspace(0, 1, pltx)
+ y = np.linspace(0, 1, pltx)
+ X, Y = np.meshgrid(x, y)
+ Z = np.zeros([pltx,pltx])
+
+ for i in range(pltx):
+ for j in range(pltx):
+ Z[j][i] = cfunc(cc,tri,K,points,[x[i],y[j]],oriented_edges,structure,edge_to_index) # compute value of pixel
+ # Z[j][i] = grad_cfunc(cc,tri,K,points,[x[i],y[j]])[0]
+
+ # plot values Z
+ fig = plt.figure()
+ plt.pcolormesh(X,Y,Z,cmap = 'jet',vmin = clim[0],vmax = clim[1])
+ plt.xlabel('x1')
+ plt.ylabel('x2')
+ plt.colorbar()
+ time = n*ht
+ title = 'test1_fineness'+str(fineness)+'_'+interpol+'_c at time step'+str(n)
+ plt.title(title)
+ image_title = title+'.png'
+ plt.savefig(image_title, bbox_inches='tight')
+ plt.close()
+ # plt.show()
+
+
+def get_linv(K_dual,nppoints_dual,numtri_dual,tri_dual,v,nodal_avrg_choice) :
+# input: triangulation (K,points,numtri etc), v=grad_c, nodal choice for linear interpolation
+# output: approximation of gradient of v
+
+ grad_cx = []
+ grad_cy = []
+ for i in range(len(K_dual)) :
+ mid = 1/3*(K_dual[i][0]+K_dual[i][1]+K_dual[i][2]) # get point inside K_dual[i]
+
+ # v is constant on K_dual, so just evalute in the middle of the triangle
+ grad_cx.append(v(mid)[0])
+ grad_cy.append(v(mid)[1])
+
+ grad_cx = np.array(grad_cx)
+ grad_cy = np.array(grad_cy)
+
+ # linearly interpolate componentwise subject to nodal_avrg_choice
+ lin_vx = getq0(K_dual,nppoints_dual,numtri_dual,tri_dual,grad_cx,nodal_avrg_choice)
+ lin_vy = getq0(K_dual,nppoints_dual,numtri_dual,tri_dual,grad_cy,nodal_avrg_choice)
+
+ return [lin_vx,lin_vy]
+
+def eval_linv(lin_vx,lin_vy,x,oriented_edges,structure,edge_to_index) :
+# input : linearly inteprolated grad_c, point x, and objects to make find_simplex for dual mesh more efficient
+# output: point evaluation of linear interpolation of piecewise constant v = grad_c
+
+ global indexE
+ indexE = find_simplex1(x,oriented_edges,structure,indexE)
+ i = edge_to_index[indexE] # get index of triangle in which x lives
+
+ val = [lin_vx[3*i]*x[0]+lin_vx[3*i+1]*x[1]+lin_vx[3*i+2], lin_vy[3*i]*x[0]+lin_vy[3*i+1]*x[1]+lin_vy[3*i+2]] # compute val=a*x+b*y+c componentwise
+
+ return val
+
+
+def approx_Laplace_cn(lin_vx,lin_vy,x,oriented_edges,structure,edge_to_index) :
+# input: normal vector, linear interpolation to piecewise constant v=grad_c, objects that make find_simplex for dual mesh more efficient
+# output: approximation to Laplace(c)
+
+ global indexE
+ indexE = find_simplex1(x,oriented_edges,structure,indexE)
+ i = edge_to_index[indexE]
+
+ Laplace = lin_vx[3*i]+lin_vy[3*i+1] # get divergence of linear interpolation of piecewise constant v=grad_c -> approximation to Laplace(c)
+
+ return Laplace
+
+
+
+# BACTERIAL DENSITY RELATED FUNCTIONS
+
+def getvKE(K,E,v) :
+# need this for functions finitevolumescheme_rho(...)
+# input: triangle K, edge E, function v
+# output: coefficient vKE from nicaise-paper, area of edge E
+
+ normalKE = unitouternormal(K,E)
+
+ aux = lambda x : np.dot(v(x),normalKE)
+
+ [integral,vKE,areaE] = avrgE(E,aux) # vKE = avrg over edge E
+
+ return [vKE,areaE]
+
+def getdistE(K1,K2) :
+# need this for function finitevolumescheme_rho(...)
+# input : neigboring triangles K1 and K2
+# output : distance of their circumcenters
+
+ CC1 = circumcenter(K1)
+ CC2 = circumcenter(K2)
+ val = np.linalg.norm(CC1-CC2, ord = None) # distance
+
+ return val
+
+def getrho(cc_old,f,ht,v,tri,K,numtri):
+# input: FE sol. c from last time step as initial value for newton method, rhs f for manufactured solutions, temporal step size ht, v=grad_c, triangulation
+# output: FV solution rho
+
+ eps = 1 # parameter epsilon = 1 in Keller-Segel system
+
+ rho = finitevolumescheme_rho(cc_old,ht,numtri,tri,K,eps,v,f)
+
+ return rho
+
+def finitevolumescheme_rho(u_old,ht,numtri,tri,K,eps,v,f) :
+# finite volume scheme to get rho with zero Neumann bdr via solving non-linear system of equations
+# input: initial value for newton method u_old, temporal step size ht, triangulation, parameter eps = 1, v=grad_c, rhs f for manufactured solutions
+# output: FV solution rho
+
+ # initialize global objects
+ global global_u_old, global_ht
+ global_u_old = u_old
+ global_ht = ht
+
+ # get non-linear system of equation to calculate num. sol. uh via FV scheme
+ def func(u) :
+
+ sum_EK = np.zeros(len(K))
+ for i in range(len(K)) : # go through triangles
+
+ [fK,avrg,areaK] = avrgK(K[i],f) # get RHS fK subject to manufactured solutions
+
+ neighbors = tri.neighbors[i] # get neighbors of triangle i
+
+ indices = [[1,2],[0,2],[0,1]] # possible combinations of vertices to compose an edge, where j = 0,1,2 is not contained , order is impotant
+
+ for j in range(3) : # go through edges
+
+ E = K[i][indices[j]] # choose edge that is shared by K[i] and K[neighbors[j]]
+
+ if neighbors[j] >= 0 : # i.e. edge E is interior
+
+ [vKE,areaE] = getvKE(K[i],E,v) # avrg over v*nKE ( = 0 if E on boundary due to neumann boundary w.r.t. c )
+ dE = getdistE(K[i],K[neighbors[j]]) # get distance between circumcenters of triangles
+
+ if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case centered flux to avoid divide by 0 or log(0)
+ sum_EK[i] = sum_EK[i] - eps*areaE/dE*(u[neighbors[j]]-u[i]) + vKE*areaE*1/2*(u[i]+u[neighbors[j]])
+ else :
+ sum_EK[i] = sum_EK[i] - eps*areaE/dE*(u[neighbors[j]]-u[i]) + vKE*areaE*((u[i]-u[neighbors[j]])/(np.log(u[i])-np.log(u[neighbors[j]])))
+
+ else : # i.e. edge E is part of the boundary -> vKE = 0 for Neumann bdr condition
+
+ sum_EK[i] = sum_EK[i] - eps*(0 + 0) #+ 0 (D_K,E) + 0 (vKE)
+
+ sum_EK[i] = 1/areaK*sum_EK[i] - 1/areaK*fK # finite volume scheme
+
+ return (u - global_u_old)/global_ht + sum_EK # difference quotient + finite volume scheme
+
+ # get Jacobian to solve non-linear system of equations for the newton method
+ def grad_func(u) :
+
+ #sum_EK = np.zeros(len(K))
+ row = []
+ col = []
+ data = []
+ for i in range(len(K)) :
+
+ [fK,avrg,areaK] = avrgK(K[i],f) # need areaK for system of equations
+ neighbors = tri.neighbors[i] # get neighbors of triangle i
+
+ indices = [[1,2],[0,2],[0,1]] # possible combinations of vertices to compose an edge, where j = 0,1,2 is not contained , order is impotant
+
+ valK = 0
+ for j in range(3) : # go through edges
+ valL = 0
+
+ E = K[i][indices[j]] # choose edge that is shared by K[i] and K[neighbors[j]]
+
+ if neighbors[j] >= 0 : # i.e. edge E is interior
+
+ [vKE,areaE] = getvKE(K[i],E,v) # avrg over v*nKE
+ dE = getdistE(K[i],K[neighbors[j]]) # get distance between circumcenters of triangles
+
+ if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case centered flux to avoid divide by 0 or log(0)
+ valK = valK + eps*(areaE/dE) + vKE*areaE*1/2
+ valL = valL - eps*(areaE/dE) + vKE*areaE*1/2
+ else :
+ valK = valK + eps*(areaE/dE) + vKE*areaE*((np.log(u[i])-np.log(u[neighbors[j]])-1+u[neighbors[j]]/u[i])/((np.log(u[i])-np.log(u[neighbors[j]]))**2))
+ valL = valL - eps*(areaE/dE) + vKE*areaE*((-np.log(u[i])+np.log(u[neighbors[j]])-1+u[i]/u[neighbors[j]])/((np.log(u[i])-np.log(u[neighbors[j]]))**2))
+
+ # get sparse Jacobian
+ row.append(i)
+ col.append(neighbors[j])
+ data.append(1/areaK*valL)
+
+
+ else : # i.e. edge E is part of the boundary -> vKE = 0 for Neumann bdr condition
+
+ valK = valK + 0
+ valL = valL + 0
+
+
+ #diagonal entries
+ row.append(i)
+ col.append(i)
+ data.append(1/global_ht+1/areaK*valK)
+
+ sparseJ = csr_matrix((data, (row, col)), shape = (numtri, numtri))#.toarray()
+
+ return sparseJ # return jacobian in sparse form
+
+ # APPLY NEWTON METHOD TO FIND ROOT U OF FUNC SUBJECT OT ITS JACOBIAN GRAD_U
+ uh = newton_method(func,grad_func, x0 = u_old ,tol = 1.5*10**-8, maxiter = 20) # tol like in fsolve
+
+ return uh
+
+
+def get_rho_plot(pltx,tri,n,ht,rho,rholim,fineness):
+# plot the finite volume solution
+
+# input: resolution of plot pltx, triangulation, time step n, step size ht, FV solution rho, plot limits rholim, fineness
+# output: plot of piecewise constatn FV sol. rho
+
+ # initialize space discretization as pixels
+ x = np.linspace(0, 1, num = pltx)
+ y = np.linspace(0, 1, num = pltx)
+ X,Y = np.meshgrid(x,y)
+ Z = np.zeros([pltx,pltx])
+
+ for i in range(pltx):
+ for j in range(pltx):
+ pt = [x[i],y[j]]
+ k = tri.find_simplex(pt) # get triangle in which pt lies
+ Z[j][i] = rho[k] # get function value
+
+ #colormap plot
+ fig = plt.figure()
+ plt.pcolormesh(X,Y,Z,cmap = 'jet',vmin = rholim[0],vmax = rholim[1]) # limits of plot
+ plt.xlabel('x1')
+ plt.ylabel('x2')
+ plt.colorbar()
+
+ # # 3D-PLOT:
+ # X, Y = np.meshgrid(x, y)
+ # fig = plt.figure()
+ # ax = plt.axes(projection='3d')
+ # plt.xlabel('x1')
+ # plt.ylabel('x2')
+ # ax.axes.set_xlim3d(left=0, right=1)
+ # ax.axes.set_ylim3d(bottom=0, top=1)
+ # ax.axes.set_zlim3d(bottom=rholim[0], top=rholim[1])
+ # p = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='jet', edgecolor='none',vmin = rholim[0], vmax = rholim[1])
+ # fig.colorbar(p, ax=ax)
+
+ time = n*ht
+ title = 'fineness'+str(fineness)+'_FV-solution_ bacterial density rho at time '+"%f" % time
+ # title_theta = 'theta_K at time t='+"%f" % time # need this for plot_aposti_estimator
+ title_png = 'FV solution rho_h at time t='+"%f" % time
+ plt.title(title_png)
+ image_title = title+'.png'
+ plt.savefig(image_title, bbox_inches='tight')
+ plt.close()
+ # plt.show()
+
+
+
+# INTERPOLATION/RECONSTRUCTION
+
+def getinterpolationRHS(tri,K,uh) :
+# reconstruct the diffusive flux of the FV scheme D_K,E = FKED (nicaise-paper notation)
+
+# input : triangulation related values tri,K and numercial solution uh
+# output : RHS for the interpolation consisting of diffusive numerical flux D_K,E = FKED (nicaise-paper notation)
+
+ FKED = [] # nicaise paper notation
+
+ for i in range(len(K)) :
+ neighbors = tri.neighbors[i] # get neighbors of triangle subject to index i
+
+ indices = [[1,2],[0,2],[0,1]] # possible combinations of vertices to compose an edge, where j = 0,1,2 is not contained , order is impotant
+
+ auxFKED = []
+ for j in range(3) : # go through neighboring triangles K[neighbors[j]] / go through edges of K[i]
+
+ E = K[i][indices[j]] # choose edge that is shared by K[i] and K[neighbors[j]]
+
+ [integral,avrg,areaE] = avrgE(E,lambda x: x) # need areaE
+
+ if neighbors[j] >= 0 : # i.e. edge E is interior
+
+ dE = getdistE(K[i],K[neighbors[j]]) # get distance between circumcenters of these two triangles
+
+ auxFKED.append(areaE/dE*(uh[neighbors[j]]-uh[i]))
+
+ else : # i.e. edge E is part of the boundary
+ auxFKED.append(0) # grad_rho*normal = 0 due to zero neumann bdr
+
+ FKED.append(auxFKED)
+
+ return FKED
+
+
+def getq0(K,nppoints,numtri,tri,uh,string: str) :
+# get linear interpolation part of the morley inteprolation q0
+
+# input : triangles K, number of triangles numtri, numerical solution uh, string on how to compute nodal values
+# output : parameters qq0 of affine linear functions q0 in form [a1,b1,c1,a2,b2,c2,...,aN,bN,cN] with N = len(K)
+
+ # initilize objects
+ rhs = []
+ row = []
+ col = []
+ data = []
+
+ for i in range(len(K)) : # go through all triangles
+
+ for j in range(3): # go through all vertices
+
+ pindex = tri.simplices[i][j] # get index of vertex
+
+ neighborstri = find_neighbors(pindex, tri) # indices of triangles that touch vertex associated with pindex
+
+ #calculate vertex_val
+ if string == 'least-squares' : # weighted mean with inner angles as weights -> see Batina, Rausch, Yang - Paper
+
+ Rx, Ry, Ixx, Iyy, Ixy = 0, 0, 0, 0, 0
+ x0 = nppoints[pindex]
+ for k in neighborstri : # go through all touching triangles
+
+ CC = circumcenter(K[k]) # calculate the cell center
+
+ Rx += (CC[0]-x0[0])
+ Ry += (CC[1]-x0[1])
+ Ixx += (CC[0]-x0[0])**2
+ Iyy += (CC[1]-x0[1])**2
+ Ixy += (CC[0]-x0[0])*(CC[1]-x0[1])
+
+ lambdax = (Ixy*Ry-Iyy*Rx)/(Ixx*Iyy-Ixy**2) # compute Lagrange multipliers
+ lambday = (Ixy*Rx-Ixx*Ry)/(Ixx*Iyy-Ixy**2)
+
+ wi = []
+ qi = []
+ for k in neighborstri : # go through all touching triangles
+ CC = circumcenter(K[k])
+ wi.append(1+lambdax*(CC[0]-x0[0])+lambday*(CC[1]-x0[1])) # get weights
+ qi.append(uh[k])
+
+ wi = np.array(wi)
+ qi = np.array(qi)
+
+ if np.sum(wi) == 0 :
+ vertex_val = np.sum(qi)/len(qi)
+ else :
+ vertex_val = np.dot(wi,qi)/np.sum(wi)
+
+
+ elif string == 'arithmetic-mean' : # calssic, arithmetic mean
+
+ vertex_val = 0
+ for k in neighborstri : # go through neighbors
+ vertex_val += uh[k]
+ vertex_val = vertex_val/len(neighborstri) # arithmetic mean
+
+ else :
+ print('Error: invalid string input')
+
+ rhs.append(vertex_val) # assign mean value over all triangles touching vertex
+
+ # get system of equaton via a*x+b*y+c, in sparse form
+ row.append(3*i+j) # for a
+ col.append(3*i+0)
+ data.append(K[i][j][0]) # triangle i, vertex j, coordinate 0 i.e. x
+ row.append(3*i+j) # for b
+ col.append(3*i+1)
+ data.append(K[i][j][1]) # triangle i, vertex j, coordinate 1 i.e. y
+ row.append(3*i+j) # for c
+ col.append(3*i+2)
+ data.append(1)
+
+ sparseQ0 = csr_matrix((data, (row, col)), shape = (3*numtri, 3*numtri))
+ qq0 = spsolve(sparseQ0, rhs, permc_spec=None, use_umfpack=True) # solve system of equations
+
+ return qq0 # coefficients piecewise linear interpolation
+
+
+def getbetaE(K,qq0,FKED) :
+# get morley coefficients betaKE
+
+# input : all tirangles K, coefficients qq0 of q0, diffusive numerical flux DKE = FKED (nicaise-paper notation)
+# output : values betaKE of form [[three-dim-array of edges of K1],[three-dim-array of edges of K2],....,[three-dim-array of edges of KN]] where N = len(K)
+
+ betaKE = [] # initilize
+
+ for i in range(len(K)) : # go through all triangles
+
+ gradq0 = lambda x : [qq0[3*i],qq0[3*i+1]] # = [a,b] if q0|K = a*x+b*y+c
+
+ indices = [[1,2],[0,2],[0,1]] # possible combinations of vertices to compose an edge, where j = 0,1,2 is not contained , order is impotant!!
+ betaE = []
+
+ for j in range(3) : # go through neighboring triangles K[neighbors[j]] i.e. go through edges of K[i]
+ E = K[i][indices[j]] # choose edge that is shared by K[i] and K[neighbors[j]]
+ normalKE = unitouternormal(K[i],E) # get normal
+
+ # compute betaKE via closed form by plugging Morley inteprolant into degree of freedom associated with betaKE
+ aux = lambda x : np.dot(gradq0(x),normalKE) # get integral over gradq0*normal i.e. normal derivative of q0
+ [integral,avrg, areaE] = avrgE(E,aux)
+ I = integral
+
+ aux = lambda x : bubbleE(K[i],j,x)*np.dot(gradbubbleK(K[i],x),normalKE) # get integral bubbleE*gradbubbleK*normal
+ [integral,avrg, areaE] = avrgE(E,aux)
+ II = integral
+
+ aux_beta = (FKED[i][j]-I)/(II) # closed form for betaKE
+
+ betaE.append(aux_beta)
+
+ betaKE.append(betaE)
+
+ return betaKE
+
+
+def morley_interpol(K,tri,qq0,betaKE,x,indexK) :
+# input : triangulation K,tri / coefficients qq0, betaKE / point x and index of the triangles in which is lives if known
+# output : function value of gradient of morley interpolant at point x
+
+ if indexK >= -0.5 : # give possibility to indicate triangle on which we work
+ i = indexK
+ else :
+ i = tri.find_simplex(x) # get triangle K[i] in which x is contained # get triangle K[i] in which x is contained
+
+ q0 = qq0[3*i]*x[0]+qq0[3*i+1]*x[1]+qq0[3*i+2] # get linear q0 on triangle K[i]
+
+ aux = 0
+
+ for j in range(3) : # sum over all edges of K[i]
+
+ aux += betaKE[i][j]*bubbleE(K[i],j,x)*bubbleK(K[i],x) # betaKE contributions to morley
+
+ return q0+aux
+
+
+def grad_morley_interpol(K,tri,qq0,betaKE,x,indexK) :
+# input : triangulation K,tri / coefficients qq0, betaKE / point x and index of the triangles in which is lives if known
+# output : function value of gradient of morley interpolant at point x
+
+ if indexK >=- 0.5 : # give possibility to indicate triangle on which we work
+ i = indexK
+ else :
+ i = tri.find_simplex(x) # get triangle K[i] in which x is contained
+
+ II = 0
+ for j in range(3) : # go through edges
+
+ II += betaKE[i][j]*(gradbubbleE(K[i],j,x)*bubbleK(K[i],x)+gradbubbleK(K[i],x)*bubbleE(K[i],j,x)) # get beta contirbutions
+
+ return np.array([qq0[3*i],qq0[3*i+1]])+II # piecewise grad_morley = grad_q0 + beta*(...) = [a,b]+beta*(...)
+
+
+def Laplace_morley(K,tri,qq0,betaKE,x,indexK) :
+# input : triangulation K,tri / coefficients qq0, betaKE / point x and index of the triangles in which is lives if known
+# output : fucntion value of the Laplacian of morley interpolant at point x
+
+ if indexK >=- 0.5 : # give possibility to indicate triangle on which we work
+ i = indexK
+ else :
+ i = tri.find_simplex(x) # get triangle K[i] in which x is contained
+
+ II = 0
+ for j in range(3) : # go through edges
+
+ aux = lambda x : [qq0[3*i],qq0[3*i+1]] # grad_q0
+ div_grad_q0 = lambda x : get_second_derivative(aux,K,tri,x,index=i) # divergence of the component-wise linear interpolation of the piecewise constant gradient grad_q0
+ II += betaKE[i][j]*(LaplacebubbleE(K[i],j,x)*bubbleK(K[i],x)+LaplacebubbleK(K[i],x)*bubbleE(K[i],j,x)+2*np.dot(gradbubbleE(K[i],j,x),gradbubbleK(K[i],x))) # product rule on beta contributions
+
+ return div_grad_q0(x)+II # piecewise Laplace_morley = "Laplace" q0 + alpha*(...)+beta*(...)
+
+
+def get_second_derivative(grad_q0,K,tri,x,index) :
+# need this for Laplacian of Morley, which is in turn needed to compute element-wise reisduals for the a posteriori error estimates
+# input : piecewise constant grad_q0 / triangulation K,tri / point x and the index of the triangle in which it lives if known
+# output: divergence of linear interpolation of piecewise constant grad_q0
+
+ if index > -0.5 : # if primal index
+ i = index
+ else :
+ i = tri.find_simplex(x)
+
+ bx = []
+ by = []
+ for j in range(3) :
+ pindex = tri.simplices[i][j] # get index of vertex
+
+ neighborstri = find_neighbors(pindex, tri) # indices of triangles that touch vertex associated with pindex
+
+ # get nodal values
+ vertex_valx = 0
+ vertex_valy = 0
+ for k in neighborstri : # go thorugh neighboring points
+ midk = 1/3*(K[k][0]+K[k][1]+K[k][2]) # middle of triangle, as function is piecewise constant it doesent matter where we evaluate
+ vertex_valx += grad_q0(midk)[0] # x-direction
+ vertex_valy += grad_q0(midk)[1] # y-direction
+ vertex_valx = vertex_valx/len(neighborstri) # arithmetic mean
+ vertex_valy = vertex_valy/len(neighborstri) # arithmetic mean
+ bx.append(vertex_valx)
+ by.append(vertex_valy)
+
+ # get linear system of equation to solve for linear inteprolation a*x+b*y+c = val
+ M = [[K[i][0][0],K[i][0][1],1],[K[i][1][0],K[i][1][1],1],[K[i][2][0],K[i][2][1],1]]
+
+ lin_vx = np.linalg.solve(M,bx)
+ lin_vy = np.linalg.solve(M,by)
+
+ div = lin_vx[0]+lin_vy[1] # get divergence of linear inteprolation
+
+ return np.array(div)
+
+
+def getmorleyplot(n,ht,K,tri,qq0,betaKE,h,rholim,fineness) :
+# input : time step n, step size ht / triangulation K,tri,points / morley coefficients qq0, betaKE / plot mesh size h, plot limits rholim, fineness
+# output : *shows/saves plot*
+
+ def exact_rho(x,t) :
+ val = 1/(1+t)*np.exp(-25*(x[0]-1/2)**2-25*(x[1]-1/2)**2)
+ return val
+
+ x = np.linspace(0, 1, int(1/h))
+ y = np.linspace(0, 1, int(1/h))
+ X, Y = np.meshgrid(x, y)
+
+ Z = np.zeros([len(X),len(Y)])
+ for i in range(len(X)) :
+ for j in range(len(Y)) :
+ pt = [x[i],y[j]]
+ Z[j][i] = morley_interpol(K,tri,qq0,betaKE,pt,indexK=-1) # plot morley interpolant
+ #Z[j][i] = exact_rho(pt,n*ht) # plot exact solution
+ #Z[j][i] = grad_morley_interpol(K,tri,v,grad_vn,qq0,betaKE,pt,indexK=-1)[1] # plot gradient of morley in one direction
+
+ # #colormap plot
+ # fig = plt.figure()
+ # plt.pcolormesh(Z,cmap = 'jet',vmin = rholim[0],vmax = rholim[1]) # limits of plot
+ # plt.xlabel('x1')
+ # plt.ylabel('x2')
+ # plt.colorbar()
+
+ # 3D-PLOT:
+ fig = plt.figure()
+ ax = plt.axes(projection='3d')
+ plt.xlabel('x1')
+ plt.ylabel('x2')
+ ax.axes.set_xlim3d(left=0, right=1)
+ ax.axes.set_ylim3d(bottom=0, top=1)
+ ax.axes.set_zlim3d(bottom=rholim[0], top=rholim[1])
+ p = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='jet', edgecolor='none',vmin = rholim[0], vmax = rholim[1])
+ fig.colorbar(p, ax=ax)
+
+ time = n*ht
+ title = 'fineness'+str(fineness)+'_morley interpolation _ bacterial density rho at time '+"%f" % time
+ title_png = 'Morley interpolation at time '+"%f" % time
+ # title = 'Exact rho at time '+"%f" % time
+ plt.title(title_png)
+ image_title = title+'.png'
+ plt.savefig(image_title, bbox_inches='tight')
+ plt.close()
+ #plt.show()
+
+
+
+# A POSTERIORI ERROR ESTIMATES
+
+def diam(K) :
+# get approximately the diameter of triangle K, need this for various residual estimates
+
+# input: triangle K
+# output: length of the longest edge, which differs the actual diameter only by a constant
+
+ indices = [[1,2],[0,2],[0,1]]
+ val_max = 0
+ for j in range(3) : # go through edges
+ E = K[indices[j]]
+ val = np.linalg.norm(E[1]-E[0]) # get length of edge
+ if val > val_max : # update val_max if val is larger
+ val_max = val
+
+ return val_max
+
+
+def get_theta_res(cc_m1,Laplace_v_m1,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index) :
+# get local residual for each triangle K
+# input: c^n-1, Laplace(c^n-1), morley^n, morley^n+1, time step size ht_0, primal and dual triangulation, objects to make find_simplex for dual mesh more efficient
+# output: the local residual part of the a posteriori error estimator
+
+ res_0 = []
+ for i in range(len(K)) :
+
+ morley_0 = lambda x: morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get morley^n
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i) # get morley^n+1
+
+ grad_morley_0 = lambda x: grad_morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get grad_morley^n
+
+ Laplace_morley_0 = lambda x: Laplace_morley(K,tri,qq0_0,betaKE_0,x,indexK=i) # get Laplace_morley^n
+
+ grad_c_m1 = lambda x : grad_cfunc(cc_m1,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) #get grad_c^n-1
+
+ # GET L2 NORM OF LOCAL RESIDUAL
+ aux_0 = lambda x : ((morley_p1(x)-morley_0(x))/ht_0 + np.dot(grad_morley_0(x),grad_c_m1(x)) + morley_0(x)*Laplace_v_m1(x) - Laplace_morley_0(x))**2
+ [integral_0,avrg,areaK] = avrgK(K[i],aux_0)
+
+ hK = diam(K[i]) # get diameter of K
+
+ res_0.append(integral_0*hK**2) # full term
+
+ return res_0
+
+
+def get_theta_diff(qq0_0,betaKE_0,K,tri) :
+# get grad_morley jump terms for each triangle K
+
+# input: morley^n and primal triangulation
+# output: grad_morley jump terms from the a posteriori error estimator
+
+ diff0 = []
+ for i in range(len(K)) : # go through triangles
+
+ indices = [[1,2],[0,2],[0,1]]
+
+ grad_morley_0 = lambda x: grad_morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get grad_morley^n
+
+ neighbors = tri.neighbors[i] # get neighbors of triangle associated with index i
+ indices = [[1,2],[0,2],[0,1]] # possible combinations of vertices to compose an edge, where j = 0,1,2 is not contained , order is impotant
+
+ val0 = 0
+ for j in range(3) : # go through edges
+
+ E = K[i][indices[j]] # choose edge that is shared by K[i] and K[neighbors[j]]
+
+ if neighbors[j] >= 0 : # i.e. edge E is interior
+
+ grad_morley_0k = lambda x: grad_morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=neighbors[j]) # get grad_morley subject to neighboring triangle
+
+ normalE = unitouternormal(K[i],E) # get normal vector
+
+ aux0 = lambda x : np.dot(grad_morley_0k(x)-grad_morley_0(x),normalE)**2 # edge L2 norm of the jump
+ [integral0,avrg,areaE] = avrgE(E,aux0)
+
+ val0 += integral0*areaE # sum over edges
+
+ else : # i.e. edge E is part of the boundary -> vKE = 0 for Neumann bdr condition
+
+ val0 += 0 # due to homogenous Neumann boundary condition
+
+ diff0.append(val0)
+
+ return diff0
+
+
+def get_theta_time(qq0_0,qq0_p1,betaKE_0,betaKE_p1,rho_m1,rho_0,rho_p1,ht_m1,ht_0,K,tri,n) :
+# get time terms for each triangle K
+
+# input: morley^n, morley^n+1, FV sol rho^n, FV sol rho^n-1, FV sol rho^n+1, step sizes ht^n-1 and ht^n, primal triangulation, time step n
+# output: time contributions to the a posteriori error estimator
+
+ if n == 0 : # slightly different situation at early times
+
+ time0 = []
+ time1 = []
+ for i in range(len(K)) :
+ morley_0 = lambda x: morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get morley^n
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i) # get morley^n+1
+
+ aux = lambda x : ((morley_p1(x)-morley_0(x))/ht_0 - (rho_p1[i]-rho_0[i])/ht_0)**2 # L2 norm to get the first time term
+ [integral,avrg,areaK] = avrgK(K[i],aux)
+
+ time0.append(integral) # first time term
+ time1.append(0) # no second time term at early times
+
+
+ else :
+ time0 = []
+ time1 = []
+ for i in range(len(K)) :
+ morley_0 = lambda x: morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get morley^n
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i) # get morley^n+1
+
+ aux = lambda x : ((morley_p1(x)-morley_0(x))/ht_0 - (rho_p1[i]-rho_0[i])/ht_0)**2 # L2 norm to get the first time term
+ [integral,avrg,areaK] = avrgK(K[i],aux)
+
+ val = areaK*((rho_p1[i]-rho_0[i])/ht_0 - (rho_0[i]-rho_m1[i])/ht_m1)**2 # # L2 norm to get the second time term
+
+ time0.append(integral) # first time term
+ time1.append(val) # second time term
+
+ return [time0,time1]
+
+
+def get_maximum_morley(rhoFV,betaKE):
+# need this for function get_theta_omega(...)
+
+# input : FV sol rhoFV and moerley coefficients betaKE
+# output: approximated Linf(Omega)-norm of morley interpoaltion
+
+ [nodal_max,beta_max] = get_max_q0_beta(rhoFV,betaKE)
+
+ return nodal_max + beta_max # sum both maximal values
+
+
+def get_max_q0_beta(rhoFV,betaKE) :
+# need this for function get_theta_omega(...)
+# idea: get maximal possible values for q0 and for beta-terms individually and then add them and use this as an upper bound
+
+# input : FV solution rhoFV and morley coefficients betaKE
+# output: maximal possible values for q0 and for beta-terms individually
+
+ nodal_max = np.max(rhoFV) # > arithmetic means to calculate nodal values for linear interpolation
+ beta_max = np.max(np.max(np.abs(betaKE),axis=0)) # > max bubbleE*bubbleK as max(bubble) = 1
+
+ return [nodal_max,beta_max]
+
+
+def get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,n) :
+# get theta_omega term from masterthesis
+
+# input: max1 and max2 that were already calculated in the previous time step, morley^n, morley^n+1, FV sol rho^n, FV sol rho^n+1
+# output: theta_omega
+
+ if n == 0 : # slightly different situation at early times
+
+ max1 = get_maximum_morley(rho_p1,betaKE_p1) # no difference
+ max2 = get_maximum_morley(np.array(rho_0)-np.array(rho_p1),np.array(betaKE_0)-np.array(betaKE_p1)) # difference
+
+ val = (max1+max2)*max2 # = theta_omega
+
+ else :
+ max3 = max1 # update terms from previous time step
+ max4 = max2
+ max1 = get_maximum_morley(rho_p1,betaKE_p1) # no difference
+ max2 = get_maximum_morley(np.array(rho_0)-np.array(rho_p1),np.array(betaKE_0)-np.array(betaKE_p1)) # difference
+
+ val = (max1+max2)*max2+max3*max4 # = theta_omega
+
+ return [val,max1,max2] # max1 and max2 will be used to calculate theta_omega for the next time step
+
+
+def get_conv_terms(v_0,qq0_p1,betaKE_p1,rho_p1,K,tri) :
+# get convection contributions for each triangle K
+
+# input: v^n=grad_c^n, morley^n+1, FV sol rho^n+1, primal triangulation
+# output: convection contributions to the a posteriori error estimator
+
+ conv_terms = []
+ for i in range(len(K)) :
+
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i) # get morley^n+1
+
+ indices = [[1,2],[0,2],[0,1]]
+ neighbors = tri.neighbors[i] # get nieghbors of triangle subject to index i
+
+ hK = diam(K[i]) # get approximated diameter of the triangle
+
+ val = 0
+ for j in range(3) : # go through all edges
+ E = K[i][indices[j]]
+
+ normalKE = unitouternormal(K[i],E) # get normal vector subject to edge E
+
+ # COMPUTE convective term F_E/areaE
+ if np.abs(rho_p1[i]-rho_p1[neighbors[j]]) < 10**-6 or rho_p1[i] < 10**-6 or rho_p1[neighbors[j]] < 10**-6 : # in this case centered flux to avoid divide by 0 or log(0)
+ F_E_over_areaE = 1/2*(rho_p1[i]+rho_p1[neighbors[j]])
+ else : # log-mean flux
+ F_E_over_areaE = (rho_p1[i]-rho_p1[neighbors[j]])/(np.log(rho_p1[i])-np.log(rho_p1[neighbors[j]]))
+
+ aux = lambda x : (np.dot(v_0(x),normalKE)*(morley_p1(x)-F_E_over_areaE))**2 # edge L2 norm
+ [integral,avrg,areaE] =avrgE(E,aux)
+
+ val += np.sqrt(integral)*hK/np.sqrt(areaE) # sum over edges
+
+ conv_terms.append(val**2)
+
+ return conv_terms
+
+
+def get_morley0_coeff(rho_0,nppoints,numtri,tri,K) :
+# get artificial morley interpolation at time zero
+# need this for extra term at early times
+
+# input: rho_0 discrete initial datum, morley^1, FV sol rho^1, primal triangulation
+# output: coefficients of artificial morley interpolation at time t=0
+
+ qq0 = getq0(K,nppoints,numtri,tri,rho_0,string = 'least-squares') # get linear interpolation q0
+
+ FKED = getinterpolationRHS(tri,K,rho_0) # get diffusive numerical fluxes DKE = FKED (nicaise-paper notation) from FV scheme
+
+ betaKE = getbetaE(K,qq0,FKED) # get coefficients betaKE
+
+ return [qq0,betaKE]
+
+
+def get_early_time_error(K,tri,v_0,qq0_0,qq0_p1,betaKE_0,betaKE_p1) :
+# get extra term at early times for each triangle K
+
+# input: primal triangulation, v^n=grad_c^n, morley^n, morley^n+1
+# output: extra term at early times for the a posteriori error estimator
+
+ res_0 = []
+ for i in range(len(K)) :
+
+ morley_0 = lambda x: morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get morley^0
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i) # get morley^1
+
+ grad_morley_0 = lambda x: grad_morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get grad_morley^0
+ grad_morley_p1 = lambda x: grad_morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i) # get grad_morley^1
+
+ aux_0 = lambda x : np.linalg.norm((morley_p1(x)-morley_0(x))*v_0(x) + (grad_morley_p1(x) - grad_morley_0(x)))**2 # L2 norm
+ [integral_0,avrg,areaK] = avrgK(K[i],aux_0)
+
+ res_0.append(integral_0)
+
+
+ return res_0
\ No newline at end of file
diff --git a/2d/plot_aposti_estimator.py b/2d/plot_aposti_estimator.py
new file mode 100644
index 0000000000000000000000000000000000000000..2be7a9a742a45858497e682024ce01a25f70cfc4
--- /dev/null
+++ b/2d/plot_aposti_estimator.py
@@ -0,0 +1,94 @@
+import myfun as my
+import numpy as np
+import time
+import pickle
+
+# ------------------------------------------------------- PLOT A POSTERIORI ERROR ESTIMATOR ---------------------------------------------------
+# This script plots the element-wise a posteriori error estimator \theta_K to see which elements would be refined when following an
+# adaptive mesh refinement strategy.
+
+tic = time.time()
+
+# SETTINGS FOR SCHEME
+nodal_avrg_choice = 'least-squares'
+method = 'wasserstein' # using log-mean convective flux
+boundary = 'Neumann'
+interpol = 'morley'
+
+test = 'blow-up'
+fineness = 3 # number of refinements of mesh before calculating numerical solution
+Nt = 8 # number of time steps used
+
+n = 6 # time step snapshot
+
+if test == 'blow-up' :
+# CONVECTION DEOMINATED REGIME i.e. BLOW-UP
+
+ def rho0(x) :
+ val = 10**(3)*np.exp(-((x[0]-0.5)**2+(x[1]-0.5)**2)/10**(-2))
+ return val
+
+ rholim = [0,5000]
+ T = 0.0045 # time interval [0,T]
+
+else :
+ print('Warning: Wrong string input for test.')
+
+# SCHEME/PROBLEM PARAMETERS
+hx = 1/2
+ht = (T)/(Nt)
+
+# GET PRIMAL MESH :
+[tri,points,nppoints,numtri] = my.initialmesh() # initialize zerost mesh for FV approximations
+K = nppoints[tri.simplices]
+
+# refine zerost mesh subject to fineness
+for i in range(fineness) :
+ hx = 1/2*hx
+ [tri,points,nppoints,numtri] = my.refinemesh(points, K, numtri)
+ K = nppoints[tri.simplices]
+
+# print some plot/solution related quantities
+print('hx',hx)
+print('ht',ht)
+pltx = 5000
+print('pltx',pltx)
+
+# initialite object
+theta_K = []
+
+if n == 0 : # slightly different situation at early times
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_early,theta_diff0,theta_time0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
+ [theta_res1,theta_diff1,theta_time1,theta_conv1,theta_conv1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # compute \theta_K for each element K
+ for i in range(len(K)) :
+ theta_K.append(theta_early[i]+theta_res1[i]+theta_diff0[i]+theta_diff1[i]+theta_time0[i]+theta_conv1[i])
+
+
+
+else : # n=/=0
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_res0,theta_diff0,theta_time0,theta_conv0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
+ [theta_res1,theta_diff1,theta_time1,theta_conv1,theta_conv1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # compute \theta_K for each element K
+ for i in range(len(theta_res1)) :
+ theta_K.append(theta_res0[i]+theta_res1[i]+theta_diff0[i]+theta_diff1[i]+theta_time0[i]+theta_time1[i]+theta_conv0[i]+theta_conv1[i])
+
+
+theta_K = np.array(theta_K)
+
+# get plot limits for error estimator
+thetalim = [0,np.max(np.max(theta_K))]
+
+# plot theta_K as piecewise constant function
+my.get_rho_plot(pltx,tri,n,ht,theta_K,thetalim,fineness)
+
\ No newline at end of file
diff --git a/2d/plot_bubble.py b/2d/plot_bubble.py
new file mode 100644
index 0000000000000000000000000000000000000000..309a4e2fbe03dd306868bac1ffa940d53423bf53
--- /dev/null
+++ b/2d/plot_bubble.py
@@ -0,0 +1,56 @@
+import myfun as my
+import numpy as np
+from matplotlib import pyplot as plt
+
+
+#PLOT BUBBLE FUNCTIONS
+
+hplt = 200
+
+K = [[0, 0],[0, 1],[1 ,0]]
+L = [[1, 1],[0, 1],[1 ,0]]
+x = np.linspace(0, 1, hplt)
+y = np.linspace(0, 1, hplt)
+X, Y = np.meshgrid(x, y)
+
+
+Z = np.zeros([len(y),len(x)])
+for i in range(len(x)) :
+ for j in range(len(y)) :
+
+ pt = [x[i],y[j]]
+ l = 0
+ if pt[0]+pt[1]<1 :
+ # Z[j][i] = my.bubbleE(K,l,pt)
+ Z[j][i] = my.bubbleK(K,pt)
+ # else:
+ # Z[j][i] = my.bubbleE(L,l,pt)
+
+
+
+# #colormap plot
+fig = plt.figure()
+plt.pcolormesh(Z,cmap = 'jet') # limits of plot
+plt.xlabel('x1')
+plt.ylabel('x2')
+plt.colorbar()
+
+# 3D-PLOT:
+# X, Y = np.meshgrid(x, y)
+# fig = plt.figure()
+# ax = plt.axes(projection='3d')
+# plt.xlabel('x1')
+# plt.ylabel('x2')
+# # ax.axes.set_xlim3d(left=0, right=1)
+# # ax.axes.set_ylim3d(bottom=0, top=1)
+# # ax.axes.set_zlim3d(bottom=rholim[0], top=rholim[1])
+# p = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='jet', edgecolor='none')
+# fig.colorbar(p, ax=ax)
+# title = 'bubbleK'
+# plt.title(title)
+# image_title = title+'.png'
+#plt.savefig(image_title, bbox_inches='tight')
+# #plt.close()
+plt.show()
+
+
\ No newline at end of file
diff --git a/2d/plot_meshes.py b/2d/plot_meshes.py
new file mode 100644
index 0000000000000000000000000000000000000000..dcceb29d72fbc7fb86824db5e875908d2d618a67
--- /dev/null
+++ b/2d/plot_meshes.py
@@ -0,0 +1,47 @@
+import myfun as my
+import numpy as np
+from scipy.spatial import Voronoi, voronoi_plot_2d
+from matplotlib import pyplot as plt
+
+
+# ------------------------------------------------------- PLOT A POSTERIORI ERROR ESTIMATOR ---------------------------------------------------
+# This script plots the primal mesh, the voronoi diagram of the primal mesh and the dual mesh in order to better understand the construction
+# of the dual mesh.
+
+
+fineness = 1 # number of refinements of mesh before calculating numerical solution
+
+[tri,points,nppoints,numtri] = my.initialmesh() # initialize zerost mesh for FV approximations
+K = nppoints[tri.simplices]
+
+hx = 1/2 # mesh size of the zerost mesh
+
+# refine zerost mesh subject to fineness
+for i in range(fineness) :
+ hx = 1/2*hx
+ [tri,points,nppoints,numtri] = my.refinemesh(points, K, numtri)
+ K = nppoints[tri.simplices]
+print('numtri',numtri)
+
+vor = Voronoi(nppoints)
+
+# get dual mesh
+[K_dual, tri_dual, points_dual, nppoints_dual, numtri_dual] = my.getdualmesh(points, tri, fineness, K)
+K_dual = np.array(K_dual)
+[tri_dual,K_dual] = my.orientation_dualmesh(tri_dual,K_dual) # enforce positive orientation of all triangles in dual mesh
+
+
+# PLOT THE PRIMAL MESH
+plt.triplot(nppoints[:,0], nppoints[:,1], tri.simplices,color='steelblue')
+plt.show()
+
+# PLOT THE VORONOI DIAGRAM, I.E. THE DIRICHLET TESSELATION, OF THE PRIMAL MESH
+fig = voronoi_plot_2d(vor)
+plt.show()
+
+# PLOT THE DUAL MESH
+plt.triplot(nppoints_dual[:,0], nppoints_dual[:,1], tri_dual.simplices,color='green')
+
+
+# plt.savefig('test.png', bbox_inches='tight')
+plt.show()
diff --git a/2d/plot_num-sol.py b/2d/plot_num-sol.py
new file mode 100644
index 0000000000000000000000000000000000000000..42c3181caca8ab80254482394fd2985811ca2e0c
--- /dev/null
+++ b/2d/plot_num-sol.py
@@ -0,0 +1,146 @@
+import myfun as my
+import numpy as np
+import time
+import pickle
+
+
+# ----------------------------------------------------------- PLOT NUMERICAL SOLUTIONS ---------------------------------------------------
+# This script plots the numerical solutions and Morley interpolations obtained via the FV-FE algorithm.
+
+tic = time.time()
+
+# SETTINGS FOR SCHEME
+nodal_avrg_choice = 'arithmetic-mean' # 'arithmetic-mean' or 'least-squares'
+method = 'wasserstein' # using log-mean convective flux
+boundary = 'Neumann'
+interpol = 'morley'
+
+test = 'manufactured1' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
+fineness = 3 # number of refinements of mesh before calculating numerical solution
+Nt = 8 # number of time steps used
+
+# initialize different test cases
+if test == 'blow-up' :
+# CONVECTION DEOMINATED REGIME i.e. BLOW-UP
+
+ def rho0(x) :
+ val = 10**(3)*np.exp(-((x[0]-0.5)**2+(x[1]-0.5)**2)/10**(-2))
+ return val
+
+ rholim = [0,5000]
+
+ T = 0.0045 # time interval [0,T]
+
+elif test == 'diff' :
+# DIFFUSION DOMINATED REGIME
+ def rho0(x) :
+ val = 1.3*np.exp(-25*(x[0]-1/2)**2-25*(x[1]-1/2)**2)
+ return val
+
+ rholim = [0,1.3]
+
+ T = 0.05 # time interval [0,T]#
+
+elif test == 'manufactured' :
+# MANUFACTURED SOLUTION WITH T = 0.5
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ rholim = [0,1.3]
+
+ T = 0.5 # time interval [0,T]#
+
+elif test == 'manufactured1' :
+# MANUFACTURED SOLUTION WITH T = 3
+
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ rholim = [0,1.3]
+
+ T = 3 # time interval [0,T]#
+
+else :
+ print('Warning: Wrong string input for test.')
+
+# SCHEME/PROBLEM PARAMETERS
+hx = 1/2
+ht = (T)/(Nt)
+
+# GET PRIMAL MESH :
+[tri,points,nppoints,numtri] = my.initialmesh() # initialize zerost mesh for FV approximations
+K = nppoints[tri.simplices]
+
+# refine zerost mesh subject to fineness
+for i in range(fineness) :
+ hx = 1/2*hx
+ [tri,points,nppoints,numtri] = my.refinemesh(points, K, numtri)
+ K = nppoints[tri.simplices]
+
+# # NEED THIS ONLY TO PLOT CHEMICAL CONCENTRATION c
+# #get dual mesh for FE method to approximate chemical density
+# [K_dual, tri_dual, points_dual, nppoints_dual, numtri_dual] = my.getdualmesh(points, tri, fineness, K)
+# K_dual = np.array(K_dual)
+# [tri_dual,K_dual] = my.orientation_dualmesh(tri_dual,K_dual) # enforce positive orientation of all triangles in dual mesh
+
+# # compute mesh topology for faster find_simplex routine for dual mesh
+# [oriented_edges,structure,edge_to_index] = my.get_edges_structure(K_dual)
+# my.init_indexE(numtri_dual)
+
+# print some plot/solution related quantities
+print('hx',hx)
+print('ht',ht)
+pltx = 200
+print('pltx',pltx)
+
+for n in [2,4,6,8] :# time snapshots that we want to plot # [0, 5, 10, 15] in kwon paper
+
+ progress = n/Nt*100
+ print("%.2f" % progress, 'procent of progress made.')
+
+ # # PLOT CHEMICAL CONCENTRATION c
+ # # load data related to c
+ # pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_c at time step'+str(n)+'.p'
+ # cc = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ # clim = [0,1.3] # set plot limits
+ # # plot piecewise linear FE solution c
+ # my.get_c_plot(pltx,tri_dual,K_dual,points_dual,n,ht,cc,clim,fineness,interpol,oriented_edges,structure,edge_to_index)
+
+ # PLOT BACTERIAL DENSITY rho
+ # load data related to rho
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_rho at time step'+str(n)+'.p'
+ [ht,rho,qq0,betaKE] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ if n == 0 :
+ [rhoFV,rho0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ # plot FV solution rho0 aka. projected initial datum for n=0
+ my.getmorleyplot(n,ht,K,tri,qq0,betaKE,1/pltx,rholim,fineness)
+ #my.get_rho_plot(pltx,tri,n,ht,rhoFV,rholim,fineness)
+ else :
+ [ht,rho,qq0,betaKE] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ # plot interpolation of Morley-type
+ my.getmorleyplot(n,ht,K,tri,qq0,betaKE,1/pltx,rholim,fineness)
+ # plot FV solution rho0
+ # my.get_rho_plot(pltx,tri,n,ht,rho,rholim,fineness)
+
+progress = 100
+print("%.2f" % progress, 'procent of progress made.')
+elapsed = time.time() - tic
+print('This took',"%.2f" % round(elapsed/60, 2), 'minutes.')
\ No newline at end of file
diff --git a/2d/scaling_space_errorestimates.py b/2d/scaling_space_errorestimates.py
new file mode 100644
index 0000000000000000000000000000000000000000..25b6a43e6b8e2610f46c6a40fc73eef768aa071e
--- /dev/null
+++ b/2d/scaling_space_errorestimates.py
@@ -0,0 +1,301 @@
+# import myfun as my
+import numpy as np
+from matplotlib import pyplot as plt
+import time
+import pickle
+
+# ----------------------------------------------------------- SCALING BEHAVIOR IN SPACE ---------------------------------------------------
+# This script takes the previously computed a posteriori error estimates on each element and adds them in the way in which they constitute
+# an upper bound for the actual error in LinfL2-L2H1 norm. Here we fix a number of time steps and investigate the scaling behavior in space,
+# i.e. making the maximal mesh size small
+
+tic = time.time()
+
+# SETTINGS FOR SCHEME
+nodal_avrg_choice = 'least-squares'
+method = 'wasserstein' # using log-mean convective flux
+boundary = 'Neumann'
+interpol = 'morley'
+
+test = 'manufactured' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
+Nt = 8 # fixed number of time steps
+
+steps = 5 # number of space refinements used
+
+# initialize different test cases
+if test == 'blow-up' :
+# CONVECTION DEOMINATED REGIME i.e. BLOW-UP
+
+ def rho0(x) :
+ val = 10**(3)*np.exp(-((x[0]-0.5)**2+(x[1]-0.5)**2)/10**(-2))
+ return val
+
+ T = 0.0045 # time interval [0,T]
+
+elif test == 'diff' :
+# DIFFUSION DOMINATED REGIME
+ def rho0(x) :
+ val = 1.3*np.exp(-25*(x[0]-1/2)**2-25*(x[1]-1/2)**2)
+ return val
+
+ T = 0.05 # time interval [0,T]#
+
+elif test == 'manufactured' :
+# MANUFACTURED SOLUTION WITH T = 0.5
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ T = 0.5 # time interval [0,T]#
+
+elif test == 'manufactured1' :
+# MANUFACTURED SOLUTION WITH T = 3
+
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ T = 3 # time interval [0,T]#
+
+else :
+ print('Warning: Wrong string input for test.')
+
+
+# SCHEME/PROBLEM PARAMETERS
+h = 1
+ht = (T)/(Nt)
+
+print('ht',ht)
+
+# initialize objects
+hh = []
+LinfL2_max = []
+L2H1_max = []
+theta_early_max = []
+theta_res_max = []
+theta_diff_max = []
+theta_time_max = []
+theta_time_val_max = []
+theta_conv_max = []
+theta_Omega_max = []
+
+
+for fineness in range(steps) : # varibale in space
+ h = 1/2*h
+ hh.append(h)
+
+ # load "exact" error of manufactured solution
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_exact error.p'
+ [L2_array,LinfL2,L2H1] = pickle.load(open('pickle_files/'+pickle_name,'rb')) # load data
+ LinfL2_max.append(LinfL2)
+ L2H1_max.append(L2H1)
+
+ # initialize objects
+ theta_early_sum = 0
+ theta_res_sum = 0
+ theta_diff_sum = 0
+ theta_time_sum = 0
+ theta_time_val_sum = 0
+ theta_conv_sum = 0
+ theta_Omega_sum = 0
+
+ for n in range(Nt) : # go through intervals inbetween time steps
+
+ if n == 0 :
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_early,theta_res0,theta_diff0,theta_time0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # initialize objects
+ theta_early_aux = 0
+ theta_res_aux = 0
+ theta_diff_aux = 0
+ theta_time_aux = 0
+ theta_time_val_aux = 0
+ theta_conv_aux = 0
+ theta_conv_aux = 0
+
+ # sum over all elements
+ for i in range(len(theta_early)) :
+ theta_early_aux += theta_early[i] # this is already the squared term
+ theta_res_aux += theta_res0[i] + theta_res1[i]
+ theta_diff_aux += theta_diff0[i]+theta_diff1[i]
+ theta_time_aux += theta_time0[i]
+ theta_conv_aux += theta_conv1[i]
+
+ # sum over time steps, corresponds to L^2 norm in time
+ theta_early_sum += ht*theta_early_aux
+ theta_res_sum += ht*theta_res_aux
+ theta_diff_sum += ht*theta_diff_aux
+ theta_time_sum += ht*theta_time_aux
+ theta_conv_sum += ht*theta_conv_aux
+ theta_Omega_sum += ht*theta_Omega0**2
+
+ elif n == Nt-1 :
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # initialize objects
+ theta_res_aux = 0
+ theta_diff_aux = 0
+ theta_time_aux = 0
+ theta_time_val_aux = 0
+ theta_conv_aux = 0
+ theta_delta_aux = 0
+
+ # sum over all elements
+ for i in range(len(theta_res1)) :
+ theta_res_aux += theta_res0[i] # this is already the squared term
+ theta_diff_aux += theta_diff0[i]
+ theta_time_aux += theta_time0[i]
+ theta_time_val_aux += theta_time_val0[i]
+ theta_conv_aux += theta_conv0[i]
+
+ # sum over time steps, corresponds to L^2 norm in time
+ theta_res_sum += ht*theta_res_aux
+ theta_diff_sum += ht*theta_diff_aux
+ theta_time_sum += ht*theta_time_aux
+ theta_time_val_sum += ht*theta_time_val_aux
+ theta_conv_sum += ht*theta_conv_aux
+ theta_Omega_sum += ht*theta_Omega0**2
+
+ else : # n=/=0
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # initialize objects
+ theta_res_aux = 0
+ theta_diff_aux = 0
+ theta_time_aux = 0
+ theta_time_val_aux = 0
+ theta_conv_aux = 0
+ theta_conv_aux = 0
+
+ # sum over all elements
+ for i in range(len(theta_res1)) :
+ theta_res_aux += theta_res0[i]+theta_res1[i] # this is already the squared term
+ theta_diff_aux += theta_diff0[i]+theta_diff1[i]
+ theta_time_aux += theta_time0[i]+theta_time1[i]
+ theta_time_val_aux += theta_time_val0[i]+theta_time_val1[i]
+ theta_conv_aux += theta_conv0[i] + theta_conv1[i]
+
+ # sum over time steps, corresponds to L^2 norm in time
+ theta_res_sum += ht*theta_res_aux
+ theta_diff_sum += ht*theta_diff_aux
+ theta_time_sum += ht*theta_time_aux
+ theta_time_val_sum += ht*theta_time_val_aux
+ theta_conv_sum += ht*theta_conv_aux
+ theta_Omega_sum += ht*theta_Omega0**2
+
+ # save results in order to compute EOC later
+ theta_early_max.append(theta_early_sum)
+ theta_res_max.append(theta_res_sum)
+ theta_diff_max.append(theta_diff_sum)
+ theta_time_max.append(theta_time_sum)
+ theta_time_val_max.append(theta_time_val_sum)
+ theta_conv_max.append(theta_conv_sum)
+ theta_Omega_max.append(theta_Omega_sum)
+
+
+# print results so far
+print('LinfL2_max ',np.round(LinfL2_max,5))
+print('L2H1_max ',np.round(L2H1_max,5))
+print('theta_early_max ',np.round(theta_early_max,5))
+print('theta_res_max ',np.round(theta_res_max,5))
+print('theta_diff_max ',np.round(theta_diff_max,40))
+print('theta_time_max ',np.round(theta_time_max,7))
+print('theta_time_val_max ',np.round(theta_time_val_max,5))
+print('theta_conv_max ',np.round(theta_conv_max,5))
+print('theta_Omega_max ',np.round(theta_Omega_max,5))
+print('mesh sizes ',hh)
+
+# initialize objects
+eoc_LinfL2 = []
+eoc_L2H1 = []
+eoc_early = []
+eoc_res = []
+eoc_diff = []
+eoc_time = []
+eoc_time_val = []
+eoc_conv = []
+eoc_Omega = []
+
+# go through all mesh sizes we used
+for i in range(steps-1) :
+
+ # COMPUTE THE EOCs
+ aux_LinfL2 = (np.log((LinfL2_max[i])/(LinfL2_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_L2H1 = (np.log((L2H1_max[i])/(L2H1_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_early = (np.log((theta_early_max[i])/(theta_early_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_res = (np.log((theta_res_max[i])/(theta_res_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_diff = (np.log((theta_diff_max[i])/(theta_diff_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_time = (np.log((theta_time_max[i])/(theta_time_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_time_val = (np.log((theta_time_val_max[i])/(theta_time_val_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_delta = (np.log((theta_conv_max[i])/(theta_conv_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_Omega = (np.log((theta_Omega_max[i])/(theta_Omega_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+
+ # save the EOCs
+ eoc_LinfL2.append(aux_LinfL2)
+ eoc_L2H1.append(aux_L2H1)
+ eoc_early.append(aux_early)
+ eoc_res.append(aux_res)
+ eoc_diff.append(aux_diff)
+ eoc_time.append(aux_time)
+ eoc_time_val.append(aux_time_val)
+ eoc_conv.append(aux_delta)
+ eoc_Omega.append(aux_Omega)
+
+
+elapsed = time.time() - tic
+print('This took',"%.2f" % round(elapsed/60, 2), 'minutes.')
+
+# print the EOCs
+print('eoc_LinfL2 ',np.round(eoc_LinfL2,5))
+print('eoc_L2H1 ',np.round(eoc_L2H1,5))
+print('eoc_early ',np.round(eoc_early,5))
+print('eoc_res ',np.round(eoc_res,6))
+print('eoc_diff ',np.round(eoc_diff,5))
+print('eoc_time ',np.round(eoc_time,5))
+print('eoc_time_val ',np.round(eoc_time_val,5))
+print('eoc_conv ',np.round(eoc_conv,5))
+print('eoc_Omega ',np.round(eoc_Omega,5))
+
+
+# plot the EOCs
+hh = np.array(hh)
+fig = plt.figure()
+plt.loglog(hh,hh,'g--',hh,hh**2,'y--',hh,LinfL2_max,'b',hh,L2H1_max,'r',hh,theta_time_max,'m',hh,theta_Omega_max,'orange',hh,theta_conv_max,'brown',hh,theta_early_max,'darkgoldenrod',hh,theta_res_max)
+
+plt.xlabel('hh')
+plt.ylabel('error estimates')
+plot_title = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error estimates'
+plt.title(plot_title)
+plt.legend(['h¹','h²','LinfL2','L2H1','theta_time','theta_Omega','theta_conv','theta_early','theta_res'])
+image_title = str(steps)+'_scaling_behavior_a-posti-error-estimates.png'#example+'_'+title+'_'+method+'_'+string+'.png'
+#plt.savefig(image_title, bbox_inches='tight')
+plt.show()
\ No newline at end of file
diff --git a/2d/scaling_time_errorestimates.py b/2d/scaling_time_errorestimates.py
new file mode 100644
index 0000000000000000000000000000000000000000..c5718be66bc2ad14682ea574a6b306f8e8d91128
--- /dev/null
+++ b/2d/scaling_time_errorestimates.py
@@ -0,0 +1,301 @@
+import numpy as np
+import myfun as my
+#from scipy.sparse.linalg import spsolve
+#from mpl_toolkits import mplot3d
+from matplotlib import pyplot as plt
+import time
+import pickle
+
+# ----------------------------------------------------------- SCALING BEHAVIOR IN SPACE ---------------------------------------------------
+# This script takes the previously computed a posteriori error estimates on each element and adds them in the way in which they constitute
+# an upper bound for the actual error in LinfL2-L2H1 norm. Here we fix a spatial mesh size and investigate the scaling behavior in time, i.e.
+# we make the temporal step size small.
+
+tic = time.time()
+
+# SETTINGS FOR SCHEME
+nodal_avrg_choice = 'least-squares'
+method = 'wasserstein' # using log-mean convective flux
+boundary = 'Neumann'
+interpol = 'morley'
+
+test = 'manufactured1' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
+
+
+# initialize different test cases
+if test == 'blow-up' :
+# CONVECTION DEOMINATED REGIME i.e. BLOW-UP
+
+ def rho0(x) :
+ val = 10**(3)*np.exp(-((x[0]-0.5)**2+(x[1]-0.5)**2)/10**(-2))
+ return val
+
+ T = 0.0045 # time interval [0,T]
+
+elif test == 'diff' :
+# DIFFUSION DOMINATED REGIME
+ def rho0(x) :
+ val = 1.3*np.exp(-25*(x[0]-1/2)**2-25*(x[1]-1/2)**2)
+ return val
+
+ T = 0.05 # time interval [0,T]#
+
+elif test == 'manufactured' :
+# MANUFACTURED SOLUTION WITH T = 0.5
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ T = 0.5 # time interval [0,T]#
+
+elif test == 'manufactured1' :
+# MANUFACTURED SOLUTION WITH T = 3
+
+ def rho0(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_rho(x,t) :
+ val = 1.3/(1+t)*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ def exact_c(x) :
+ val = 1.3*np.exp(-60*(x[0]-1/2)**2-60*(x[1]-1/2)**2)
+ return val
+
+ T = 3 # time interval [0,T]#
+
+else :
+ print('Warning: Wrong string input for test.')
+
+
+# SCHEME/PROBLEM PARAMETERS
+
+# initialize objects
+hh = []
+LinfL2_max = []
+L2H1_max = []
+theta_early_max = []
+theta_res_max = []
+theta_diff_max = []
+theta_time_max = []
+theta_time_val_max = []
+theta_conv_max = []
+theta_delta_max = []
+theta_Omega_max = []
+
+fineness = 4 # fix a number of refinements to the zerost mesh, i.e. fix a spatial mesh size
+steps = 8 # number of different time step sizes used
+
+for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4, 6 ,8 , 16, 32, 64]
+
+ ht = T/Nt
+ hh.append(ht)
+
+ # initialize objects
+ theta_early_sum = 0
+ theta_res_sum = 0
+ theta_diff_sum = 0
+ theta_time_sum = 0
+
+ theta_time_val_sum = 0
+ theta_conv_sum = 0
+ theta_delta_sum = 0
+ theta_Omega_sum = 0
+
+ # load 'exact' error
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_exact error.p'
+ [L2_array,LinfL2,L2H1] = pickle.load(open('pickle_files/'+pickle_name,'rb')) # load data
+ LinfL2_max.append(LinfL2)
+ L2H1_max.append(L2H1)
+
+ for n in range(Nt) : # go through intervals inbetween time steps
+
+ if n == 0 :
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_early,theta_res0,theta_diff0,theta_time0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # initialize objects
+ theta_early_aux = 0
+ theta_res_aux = 0
+ theta_diff_aux = 0
+ theta_time_aux = 0
+ theta_conv_aux = 0
+ theta_delta_aux = 0
+
+ # sum over all elements
+ for i in range(len(theta_early)) :
+ theta_early_aux += theta_early[i] # this is already the squared term
+ theta_res_aux += theta_res0[i] + theta_res1[i]
+ theta_diff_aux += theta_diff0[i]+theta_diff1[i]
+ theta_time_aux += theta_time0[i]
+ theta_conv_aux += theta_conv1[i]
+
+ # sum over time steps, corresponds to L^2 norm in time
+ theta_early_sum += ht*theta_early_aux
+ theta_res_sum += ht*theta_res_aux
+ theta_diff_sum += ht*theta_diff_aux
+ theta_time_sum += ht*theta_time_aux
+ theta_delta_sum += ht*theta_delta_aux
+ theta_Omega_sum += ht*theta_Omega0**2
+
+ elif n == Nt-1 :
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # initialize objects
+ theta_res_aux = 0
+ theta_diff_aux = 0
+ theta_time_aux = 0
+ theta_time_val_aux = 0
+ theta_conv_aux = 0
+ theta_delta_aux = 0
+
+ # sum over all elements
+ for i in range(len(theta_res1)) :
+ theta_res_aux += theta_res0[i] # this is already the squared term
+ theta_diff_aux += theta_diff0[i]
+ theta_time_aux += theta_time0[i]
+ theta_time_val_aux += theta_time_val0[i]
+ theta_conv_aux += theta_conv0[i]
+
+ # sum over time steps, corresponds to L^2 norm in time
+ theta_res_sum += ht*theta_res_aux
+ theta_diff_sum += ht*theta_diff_aux
+ theta_time_sum += ht*theta_time_aux
+ theta_time_val_sum += ht*theta_time_val_aux
+ theta_conv_sum += ht*theta_conv_aux
+ theta_Omega_sum += ht*theta_Omega0**2
+
+ else : # n=/=0
+
+ # load a posteriori error estimates
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+
+ # initialize objects
+ theta_res_aux = 0
+ theta_diff_aux = 0
+ theta_time_aux = 0
+ theta_time_val_aux = 0
+ theta_conv_aux = 0
+ theta_delta_aux = 0
+
+ # sum over all elements
+ for i in range(len(theta_res1)) :
+ theta_res_aux += theta_res0[i]+theta_res1[i] # this is already the squared term
+ theta_diff_aux += theta_diff0[i]+theta_diff1[i]
+ theta_time_aux += theta_time0[i]+theta_time1[i]
+ theta_time_val_aux += theta_time_val0[i]+theta_time_val1[i]
+ theta_conv_aux += theta_conv0[i]+theta_conv1[i]
+
+ # sum over time steps, corresponds to L^2 norm in time
+ theta_res_sum += ht*theta_res_aux
+ theta_diff_sum += ht*theta_diff_aux
+ theta_time_sum += ht*theta_time_aux
+ theta_time_val_sum += ht*theta_time_val_aux
+ theta_conv_sum += ht*theta_conv_aux
+ theta_Omega_sum += ht*theta_Omega0**2
+
+ # save results in order to compute EOC later
+ theta_early_max.append(theta_early_sum)
+ theta_res_max.append(theta_res_sum)
+ theta_diff_max.append(theta_diff_sum)
+ theta_time_max.append(theta_time_sum)
+ theta_time_val_max.append(theta_time_val_sum)
+ theta_conv_max.append(theta_conv_sum)
+ theta_Omega_max.append(theta_Omega_sum)
+
+# print results so far
+print('LinfL2_max ',np.round(LinfL2_max,5))
+print('L2H1_max ',np.round(L2H1_max,5))
+print('theta_early_max ',np.round(theta_early_max,5))
+print('theta_res_max ',np.round(theta_res_max,5))
+print('theta_diff_max ',np.round(theta_diff_max,40))
+print('theta_time_max ',np.round(theta_time_max,6))
+print('theta_time_val_max ',np.round(theta_time_val_max,5))
+print('theta_conv_max ',np.round(theta_conv_max,5))
+print('theta_Omega_max ',np.round(theta_Omega_max,5))
+print('mesh sizes ',hh)
+
+# initialize objects
+eoc_LinfL2 = []
+eoc_L2H1 = []
+eoc_early = []
+eoc_res = []
+eoc_diff = []
+eoc_time = []
+eoc_time_val = []
+eoc_conv = []
+eoc_Omega = []
+
+# go through all numbers of number of step sizes we used
+for i in range(steps-1) :
+
+ # COMPUTE EOCs
+ aux_LinfL2 = (np.log((LinfL2_max[i])/(LinfL2_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_L2H1 = (np.log((L2H1_max[i])/(L2H1_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_early = (np.log((theta_early_max[i])/(theta_early_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_res = (np.log((theta_res_max[i])/(theta_res_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_diff = (np.log((theta_diff_max[i])/(theta_diff_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_time = (np.log((theta_time_max[i])/(theta_time_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_time_val = (np.log((theta_time_val_max[i])/(theta_time_val_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_conv = (np.log((theta_conv_max[i])/(theta_conv_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_Omega = (np.log((theta_Omega_max[i])/(theta_Omega_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+
+ # save the EOCs
+ eoc_LinfL2.append(aux_LinfL2)
+ eoc_L2H1.append(aux_L2H1)
+ eoc_early.append(aux_early)
+ eoc_res.append(aux_res)
+ eoc_diff.append(aux_diff)
+ eoc_time.append(aux_time)
+ eoc_time_val.append(aux_time_val)
+ eoc_conv.append(aux_conv)
+ eoc_Omega.append(aux_Omega)
+
+
+elapsed = time.time() - tic
+print('This took',"%.2f" % round(elapsed/60, 2), 'minutes.')
+
+# print the EOCs
+print('eoc_LinfL2 ',np.round(eoc_LinfL2,5))
+print('eoc_L2H1 ',np.round(eoc_L2H1,5))
+print('eoc_early ',np.round(eoc_early,5))
+print('eoc_res ',np.round(eoc_res,5))
+print('eoc_diff ',np.round(eoc_diff,5))
+print('eoc_time ',np.round(eoc_time,5))
+print('eoc_time_val ',np.round(eoc_time_val,5))
+print('eoc_conv ',np.round(eoc_conv,5))
+print('eoc_Omega ',np.round(eoc_Omega,5))
+
+
+# plot the EOCs
+hh = np.array(hh)
+fig = plt.figure()
+plt.loglog(hh,hh,'g--',hh,hh**2,'y--',hh,LinfL2_max,'b',hh,L2H1_max,'r',hh,theta_time_max,'m',hh,theta_Omega_max,'orange',hh,theta_conv_max,'brown',hh,theta_early_max,'darkgoldenrod',hh,theta_res_max)
+
+plt.xlabel('hh')
+plt.ylabel('error estimates')
+plot_title = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error estimates'
+plt.title(plot_title)
+plt.legend(['h¹','h²','LinfL2','L2H1','theta_time','theta_Omega','theta_conv','theta_early','theta_res'])
+image_title = str(steps)+'_scaling_behavior_a-posti-error-estimates.png'#example+'_'+title+'_'+method+'_'+string+'.png'
+#plt.savefig(image_title, bbox_inches='tight')
+# plt.show()