Notebook

This file is part of https://github.com/diehlpk/reusommer21.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.ticker import FormatStrFormatter
import sys

In [2]:

example = "Quartic"
has_condition = True
save_results = True
con = []

In [3]:

#############################################################################
# Exact solution 
#############################################################################

def exactSolution(x):
    
    if example == "Cubic":
        return (2/3/np.sqrt(3)) * ( 9*x - 9*x*x + 2 * x * x * x )
    elif example == "Quartic":
        return 16/9 * x * x - 32/27 * x * x * x + 16/81 * x * x * x * x
    elif example == "Quadratic":
        return  4/3 * x - 4/9 * x * x
    else:
        print("Error: Either provide Linear, Quadratic, Quartic, or Cubic")
        sys.exit()

In [4]:

#############################################################################
# Solve the system
#############################################################################

def solve(M,f):
    return np.linalg.solve(M,f)

In [5]:

#############################################################################
# Loading
#############################################################################

def f(x):
    
    if example == "Cubic":
        return -( 2/np.sqrt(3)) * ( -6 + 4*x )
    elif example == "Quartic":
        return  -32/9 + 64/9 * x - 64/27 * x * x
    elif example == "Quadratic":
        return 8/9
    else:
        print("Error: Either provide Quadratic, Quartic, or Cubic")
        sys.exit()

def forceFull(n,x):
    
    force = np.zeros(n)
   
    for i in range(1,n-1):
        force[i] = f(x[i])
    
    force[n-1] = 0
    
    return force

def forceCoupling(n,x):
    
    dim = 2*n + 2*n-1
    
    force = np.zeros(dim)
   
    for i in range(1,dim-1):
        force[i] = f(x[i])
    
    force[dim-1] = 0
    
    return force

In [6]:

#############################################################################
# Assemble the stiffness matrix for the finite difference model (FD)
#############################################################################

def FDM(n,h):

    M = np.zeros([n,n])

    M[0][0] = 1

    for i in range(1,n-1):
        M[i][i-1] = -2
        M[i][i] = 4
        M[i][i+1] = -2

    M[n-1][n-1] = 1

    M *= 1./(2.*h*h)

    return M

In [7]:

#############################################################################
# Assemble the stiffness matrix for the coupling of FDM - VHM - FDM
#############################################################################

def CouplingFDVHM(n,h):

    fVHM = 1./(8.*h/2*h/2)
    fFDM = 1./(2.*h*h)
    
    dim = 2*n + 2*n-1
    

    M = np.zeros([dim,dim])
    
    M[0][0] = 1 

    for i in range(1,n-1):
        M[i][i-1] = -2 * fFDM
        M[i][i] = 4 * fFDM
        M[i][i+1] = -2 * fFDM

    M[n-1][n-1] = -1 
    M[n-1][n] = 1  

    M[n][n-1] = 11*h * fFDM / 3
    M[n][n-2] = -18*h * fFDM / 3
    M[n][n-3] = 9*h * fFDM / 3
    M[n][n-4] = -2*h * fFDM / 3

    M[n][n] = 11 / 6 / h *2
    M[n][n+1] = -18 / 6 / h *2
    M[n][n+2] = 9 / 6 / h *2
    M[n][n+3] = -2 / 6 / h *2

    M[n+1][n] = -8 * fVHM
    M[n+1][n+1] = 16 * fVHM
    M[n+1][n+2] = -8 * fVHM
    
    mid = n+2*n -2

    for i in range(n+2,mid-1):
        M[i][i-2] = -1. * fVHM
        M[i][i-1] = -4. * fVHM
        M[i][i] = 10. * fVHM
        M[i][i+1] =  -4. * fVHM
        M[i][i+2] = -1. * fVHM
        
   
        
    M[mid-1][mid-2] = -8 * fVHM
    M[mid-1][mid-1] = 16 * fVHM
    M[mid-1][mid] = -8 * fVHM

    M[mid][mid+1] = -1 
    M[mid][mid] = 1  

    M[mid+1][mid+1-1] = 11 / 6 / h *2
    M[mid+1][mid+1-2] = -18 / 6 / h*2
    M[mid+1][mid+1-3] = 9 / 6 / h*2
    M[mid+1][mid+1-4] = -2 / 6 / h*2

    M[mid+1][mid+1] = 11 *h * fFDM / 3
    M[mid+1][mid+1+1] = -18 *h * fFDM / 3
    M[mid+1][mid+1+2] = 9 *h * fFDM / 3
    M[mid+1][mid+1+3] = -2 *h * fFDM / 3

    for i in range(mid+2,dim-1):
        M[i][i-1] = -2 * fFDM
        M[i][i] = 4 * fFDM
        M[i][i+1] = -2 * fFDM

    M[dim-1][dim-1] = 1
    
    if has_condition:
        con.append(np.linalg.cond(M))

    return M

In [8]:

def compute(i):
    
    n = np.power(2,int(i))
    h = 1./n
    nodes = n + 1
    
    print(nodes,h)
    x1 = np.linspace(0,1,nodes)
    x2 = np.linspace(1,2.,2*nodes-1)
    x3 = np.linspace(2,3.,nodes)
    x = np.array(np.concatenate((x1,x2,x3)))
    
    M = CouplingFDVHM(nodes,h)
    
    f = forceCoupling(nodes,x)
    
    f[n] = 0
    f[n+1] = 0
    
    mid = n+2*n+1

    f[mid] = 0
    f[mid+1] = 0
    
    u = solve(M,f)
    
    x1 = x[0:nodes-1]
    x2 = x[nodes:mid+1:2]
    x3 = x[mid+2:len(x)]
    
    u1 = u[0:nodes-1]
    u2 = u[nodes:mid+1:2]
    u3 = u[mid+2:len(x)]
    
    x = np.concatenate([x1,x2,x3])
    u = np.concatenate([u1,u2,u3])
    
    plt.plot(x,u,label="$hFD=$"+str(h))
    plt.grid(True)
    
    x = np.arange(0,3+h,h)
    interval = int(len(x)/3)
    
    forceFD = forceFull(len(x),x)
    MFD = FDM(len(x),h)
    
    uFM = solve(MFD,forceFD)
      
    uFD = uFM
    
    if save_results:
        np.savetxt("u_vhcm_d_"+str(n)+".csv", u, delimiter=",")

    return u , uFD , x

In [9]:

uCoupled = []
uLocal = []
xLocal = []
max_iteration = 8

for i in range(4,max_iteration):
    u , uFD , x = compute(i)
    uCoupled.append(u)
    uLocal.append(uFD)
    xLocal.append(x)

n = np.power(2,int(max_iteration))
h = 1./n

x = np.arange(0,3+h,h)
plt.plot(x,exactSolution(x),label="Exact solution")
plt.legend()
plt.xlabel(r"Position $x$")
plt.ylabel(r"Displacement $u$")
plt.title(str(example)+" solution ")
plt.savefig("coupling-vhm-"+str(example)+"-dirchelt.pdf")

17 0.0625
33 0.03125
65 0.015625
129 0.0078125

In [10]:

plt.clf()
markers = ['s','o','x','.']
level = [8,16,32,64]

for i in range(0,len(xLocal)):
    n = np.power(2,i+4)
    plt.plot(xLocal[i],uLocal[i]-uCoupled[i],color="black",marker=markers[i],markevery=level[i],label=r"$\delta$=1/"+str(int(n/2)))
    print('{0:10f}'.format(max(abs(uLocal[i]-uCoupled[i]))))

plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%0.6f'))
plt.xlabel(r"$x$")
plt.ylabel(r"Error in displacement w.r.t FDM")
plt.title(r"Example with " + str(example).lower() + " solution for VHCM with m = 2")
plt.grid()
plt.axvline(x=1,c="#536872")
plt.axvline(x=2,c="#536872")
plt.legend()
plt.savefig("vhm-"+str(example)+"-d-error.pdf",bbox_inches='tight')

if has_condition :
    np.savetxt("con_vhm_d.csv", con, delimiter=",")

In [ ]: