From 75981c47d1f5d5d045205fc896aa8caf586e260b Mon Sep 17 00:00:00 2001 From: Lars Gebraad Date: Thu, 22 Mar 2018 14:59:48 +0100 Subject: [PATCH] Reworked the main method of animation (now a sequence of images). Added standard configuration for time-stepping. --- alt_main.py | 113 ---------------------------------- config.py | 2 + main.py | 171 +++++++++++++++++++++++++++++++++++++++++++++++++++ operators.py | 82 +----------------------- 4 files changed, 174 insertions(+), 194 deletions(-) delete mode 100644 alt_main.py create mode 100644 main.py diff --git a/alt_main.py b/alt_main.py deleted file mode 100644 index 64b868f..0000000 --- a/alt_main.py +++ /dev/null @@ -1,113 +0,0 @@ -#!/usr/bin/python3.6 -import matplotlib.pyplot as plt -import numpy as np -import scipy.sparse as sparse -from scipy.sparse.linalg import inv -import discretizationUtils as grid -import config as config -import operators as op -import matrices -import matplotlib.animation as animation - -# --- Setup; using the config as a 'global variable module' --- # -config.dt = 0.0001 -config.ny = 50 -config.nx = 100 -config.dy = 1 / (config.ny - 1) -config.dx = 1 / (config.ny - 1) - -dt = config.dt -ny = config.ny -nx = config.nx -dy = config.dy -dx = config.dx - -# Physical regime -sqrtRa = 7 -Tbottom = 1 -Ttop = 0 -tempDirichlet = [Tbottom, Ttop] -tempNeumann = [0, 0] -# --- end --- # - -# Initial conditions -startT = np.expand_dims(np.linspace(Tbottom, Ttop, ny), 1).repeat(nx, axis=1) -startT[1:4, 9:11] = 1.1 -startT[1:4, 29:31] = 1.1 -startT = grid.fieldToVec(startT) -startPsi = np.expand_dims(np.zeros(ny), 1).repeat(nx, axis=1) # start with zeros for streamfunctions -startPsi = grid.fieldToVec(startPsi) - -# Generate operators for differentials (generating them in functions is more computationally expensive) -dyOpPsi = op.dyOp() -dxOpPsi = op.dxOpStream() -dlOpPsi = op.dlOpStreamMod() - -dlOpTemp, rhsDlOpTemp = op.dlOpTemp(bcDirArray=tempDirichlet, bcNeuArray=tempNeumann) -dxOpTemp, rhsDxOpTemp = op.dxOpTemp(bcNeuArray=tempNeumann) -dyOpTemp, rhsDyOpTemp = op.dyOpTemp(bcDirArray=tempDirichlet) - -Tnplus1 = startT -Psinplus1 = startPsi - -toKeep = np.ones((nx * ny,)) -toKeep[0:ny] = 0 -toKeep[-ny:] = 0 -toKeep[ny::ny] = 0 -toKeep[ny - 1::ny] = 0 -toKeep = sparse.csr_matrix(sparse.diags(toKeep, 0)) - -fig = plt.figure(figsize=(10, 5), dpi=400) - -Writer = animation.writers['ffmpeg'] -writer = Writer(fps=30, metadata=dict(artist='Lars Gebraad'), bitrate=5000) -ims = [] - -invDlOpPsi = sparse.linalg.inv(dlOpPsi) -print("Done!") - -plt.imshow(invDlOpPsi.todense()) -plt.show() - -# Integrate in time -for it in range(500): - Tn = Tnplus1 - Psin = Psinplus1 - C, rhsC = matrices.constructC(dxOpTemp=dxOpTemp, rhsDxOpTemp=rhsDxOpTemp, dyOpTemp=dyOpTemp, - rhsDyOpTemp=rhsDyOpTemp, - dlOpTemp=dlOpTemp, rhsDlOpTemp=rhsDlOpTemp, psi=Psin, dxOpPsi=dxOpPsi, - dyOpPsi=dyOpPsi, sqrtRa=sqrtRa) - # Tnplus1 = inv(sparse.eye(nx*ny) + (dt/2) * C) @ (dt * rhsC + (sparse.eye(nx*ny) - (dt/2) * C) @ Tn) - Tnplus1 = sparse.linalg.spsolve(sparse.eye(nx * ny) + (dt / 2) * C, - dt * rhsC + (sparse.eye(nx * ny) - (dt / 2) * C) @ Tn, use_umfpack=True) - - # Enforcing column vector - Tnplus1.shape = (nx * ny, 1) - - # Enforcing Dirichlet - Tnplus1[0::ny] = Tbottom * 2 - Tnplus1[ny - 1::ny] = Ttop - - Psinplus1 = invDlOpPsi @ (- sqrtRa * toKeep @ (dxOpTemp @ Tnplus1 - rhsDxOpTemp)) - # Psinplus1 = sparse.linalg.spsolve(dlOpPsi, - sqrtRa * toKeep @ (dxOpTemp @ Tnplus1 - rhsDxOpTemp),use_umfpack=True) - - # Enforcing column vector - Psinplus1.shape = (nx * ny, 1) - # Reset Dirichlet - Psinplus1[0::ny] = 0 - Psinplus1[ny - 1::ny] = 0 - Psinplus1[0:ny] = 0 - Psinplus1[-ny:] = 0 - - # if (it % 100 == 0): - # im = plt.imshow(np.flipud(grid.vecToField(Tnplus1)), animated=True,vmin=-1, vmax=2, aspect=1, cmap=plt.get_cmap('magma')) - # ims.append([im]) - -# plt.imshow(grid.vecToField(dyOpTemp @ startT - rhsDyOpTemp), vmin=-1, vmax=1) -plt.imshow(grid.vecToField(Tnplus1), vmin=0) -plt.colorbar() -plt.show() - -# ani = animation.ArtistAnimation(fig, ims, interval=10, blit=True, repeat_delay=0) -# # plt.colorbar() -# ani.save('animation.mp4', writer=writer) diff --git a/config.py b/config.py index 9195c4f..82a285f 100644 --- a/config.py +++ b/config.py @@ -3,3 +3,5 @@ ny = 10 dx = 0.1 dy = 0.1 +dt = 0.0001 +nt = 500 \ No newline at end of file diff --git a/main.py b/main.py new file mode 100644 index 0000000..d693dd8 --- /dev/null +++ b/main.py @@ -0,0 +1,171 @@ +#!/usr/bin/python3.6 +import matplotlib.pyplot as plt +import numpy as np +import scipy.sparse as sparse +import scipy.sparse.linalg +import discretizationUtils as grid +import config as config +import operators as op +import matrices +import matplotlib.animation as animation +import time + +# --- Setup; using the config as a 'global variable module' --- # +config.dt = 0.001 +config.nt = 1500 +config.ny = 20 +config.nx = 30 +config.dy = 1 / (config.ny - 1) +config.dx = 1 / (config.ny - 1) + +dt = config.dt +nt = config.nt +ny = config.ny +nx = config.nx +dy = config.dy +dx = config.dx + +# Physical regime +sqrtRa = 7 +Tbottom = 1 +Ttop = 0 +tempDirichlet = [Tbottom, Ttop] +tempNeumann = [0, 0] + +# Takes a long time! +generateAnimation = True +generateEvery = 20 +# Speeds up! +useSuperLUFactorizationEq1 = True +useSuperLUFactorizationEq2 = True + +# --- No modifications below this point! --- # + +# Initial conditions +startT = np.expand_dims(np.linspace(Tbottom, Ttop, ny), 1).repeat(nx, axis=1) +# startT[1:4, 9:11] = 1.1 +# startT[1:4, 29:31] = 1.1 +startT = grid.fieldToVec(startT) +startPsi = np.expand_dims(np.zeros(ny), 1).repeat(nx, axis=1) # start with zeros for streamfunctions +startPsi = grid.fieldToVec(startPsi) + +# Generate operators for differentials +dyOpPsi = op.dyOp() +dxOpPsi = op.dxOpStream() +dlOpPsi = op.dlOpStreamMod() +dlOpTemp, rhsDlOpTemp = op.dlOpTemp(bcDirArray=tempDirichlet, bcNeuArray=tempNeumann) +dxOpTemp, rhsDxOpTemp = op.dxOpTemp(bcNeuArray=tempNeumann) +dyOpTemp, rhsDyOpTemp = op.dyOpTemp(bcDirArray=tempDirichlet) + +# This matrix is needed to only use the Non-Dirichlet rows in the del²psi operator. The other rows are basically +# ensuring psi_i,j = 0. What it does is it removes rows from a sparse matrix (unit matrix with some missing elements). +psiNonDirichlet = np.ones((nx * ny,)) +psiNonDirichlet[0:ny] = 0 +psiNonDirichlet[-ny:] = 0 +psiNonDirichlet[ny::ny] = 0 +psiNonDirichlet[ny - 1::ny] = 0 +psiNonDirichlet = sparse.csc_matrix(sparse.diags(psiNonDirichlet, 0)) + +# Pretty straightforwad; set up the animation stuff +if (generateAnimation): + # Set up for animation + fig = plt.figure(figsize=(8.5, 5), dpi=300) + Writer = animation.writers['ffmpeg'] + writer = Writer(fps=30, metadata=dict(artist='Lars Gebraad'), bitrate=5000) + ims = [] +else: + fig = plt.figure(figsize=(8.5, 5), dpi=600) + +# -- Preconditioning -- # +if (useSuperLUFactorizationEq1): + start = time.time() + factor1 = sparse.linalg.factorized(dlOpPsi) # Makes LU decomposition. + end = time.time() + print("Runtime of LU factorization for eq-1: ", end - start, "seconds") + +# Set up initial +Tnplus1 = startT +Psinplus1 = startPsi + +start = time.time() +# Integrate in time +print("Starting time marching for", nt, "steps ...") +t = [0] +for it in range(nt): + # if(it == 20): + # dt = dt*10 + t.append(t[-1] + dt) + Tn = Tnplus1 + Psin = Psinplus1 + + # Regenerate C + C, rhsC = matrices.constructC(dxOpTemp=dxOpTemp, rhsDxOpTemp=rhsDxOpTemp, dyOpTemp=dyOpTemp, + rhsDyOpTemp=rhsDyOpTemp, dlOpTemp=dlOpTemp, rhsDlOpTemp=rhsDlOpTemp, psi=Psin, + dxOpPsi=dxOpPsi, dyOpPsi=dyOpPsi, sqrtRa=sqrtRa) + + # Solve for Tn+1 + if (useSuperLUFactorizationEq2): + factor2 = sparse.linalg.factorized( + sparse.csc_matrix(sparse.eye(nx * ny) + (dt / 2) * C)) # Makes LU decomposition. + Tnplus1 = factor2(dt * rhsC + (sparse.eye(nx * ny) - (dt / 2) * C) @ Tn) + else: + Tnplus1 = sparse.linalg.spsolve(sparse.eye(nx * ny) + (dt / 2) * C, + dt * rhsC + (sparse.eye(nx * ny) - (dt / 2) * C) @ Tn) + + # Enforcing column vector + Tnplus1.shape = (nx * ny, 1) + + # Enforcing Dirichlet + Tnplus1[0::ny] = Tbottom + Tnplus1[ny - 1::ny] = Ttop + + if (useSuperLUFactorizationEq1): + Psinplus1 = factor1(- sqrtRa * psiNonDirichlet @ (dxOpTemp @ Tnplus1 - rhsDxOpTemp)) + else: + Psinplus1 = sparse.linalg.spsolve(dlOpPsi, - sqrtRa * psiNonDirichlet @ (dxOpTemp @ Tnplus1 - rhsDxOpTemp)) + + # Enforcing column vector + Psinplus1.shape = (nx * ny, 1) + # Reset Dirichlet + Psinplus1[0::ny] = 0 + Psinplus1[ny - 1::ny] = 0 + Psinplus1[0:ny] = 0 + Psinplus1[-ny:] = 0 + + if (it % generateEvery == 0): + if (generateAnimation): + ux = grid.vecToField(dyOpPsi @ Psinplus1) + uy = grid.vecToField(dxOpPsi @ Psinplus1) + maxSpeed = np.max(np.square(ux) + np.square(uy)) + plt.quiver(np.flipud(ux), -np.flipud(uy)) + im = plt.imshow(np.flipud(grid.vecToField(Tnplus1)), vmin=0, aspect=1, cmap=plt.get_cmap('magma'), vmax=2) + clrbr = plt.colorbar() + plt.tight_layout() + plt.title("t: %.2f, it: %i, Maximum speed: %.2f" % (t[-1], it, maxSpeed)) + plt.savefig("fields/field%i.png" % (it / generateEvery)) + plt.gca().clear() + clrbr.remove() + # ims.append([im]) + print("time step:", it, end='\r') + +end = time.time() +print("Runtime of time-marching: ", end - start, "seconds") + +# And output figure(s) +if (generateAnimation): + # ani = animation.ArtistAnimation(fig, ims, interval=10, blit=True, repeat_delay=0) + # plt.colorbar() + # plt.tight_layout() + # plt.title("Maximum speed: %.2f" % maxSpeed) + # ani.save('animation.mp4', writer=writer) + a = 1 +else: + ux = grid.vecToField(dyOpPsi @ Psinplus1) + uy = grid.vecToField(dxOpPsi @ Psinplus1) + maxSpeed = np.max(np.square(ux) + np.square(uy)) + plt.quiver(np.flipud(ux), -np.flipud(uy)) + plt.imshow(np.flipud(grid.vecToField(Tnplus1)), cmap=plt.get_cmap('magma')) + plt.colorbar() + plt.tight_layout() + plt.title("Maximum speed: %.2f" % maxSpeed) + plt.savefig("Convection-field.png") diff --git a/operators.py b/operators.py index 43878f7..53d8761 100644 --- a/operators.py +++ b/operators.py @@ -197,86 +197,6 @@ def dlOpTemp(bcDirArray, bcNeuArray): return A, rhs -# def dlOpTempLegacy(bcDirArray, bcNeuArray): -# # TODO -# # Streamline Laplace operator -# nx = config.nx -# ny = config.ny -# dx = config.dx -# dy = config.dy -# -# rhs = np.zeros((nx * ny,)) -# rhs.shape = (nx * ny, 1) -# -# mid = -2 * np.ones((ny,)) / (dy ** 2) - 2 * np.ones((ny,)) / (dx ** 2) -# mid[0] = 0 -# mid[-1] = 0 -# mid = np.tile(mid, (nx,)) -# # These are the one-sided differences for d²T/dx² -# mid[1:ny - 1] = -2 / (dy ** 2) - 1 / (2 * dx ** 2) -# mid[-ny + 1:-1] = -2 / (dy ** 2) - 1 / (2 * dx ** 2) -# -# plus1x = np.ones((ny,)) / (dx ** 2) -# plus1x[0] = 0 -# plus1x[-1] = 0 -# plus1x = np.tile(plus1x, (nx - 1,)) -# plus1x[0:ny] = 0 -# -# min1x = np.ones((ny,)) / (dx ** 2) -# min1x[0] = 0 -# min1x[-1] = 0 -# min1x = np.tile(min1x, (nx - 1,)) -# min1x[-ny:] = 0 -# # Associated RHS -# rhs[ny::ny] = -bcDirArray[0] / (dx ** 2) -# rhs[ny +::ny] = -bcDirArray[0] / (dx ** 2) -# -# # Check -# plus2x = np.zeros((ny,)) -# plus2x = np.tile(plus2x, (nx - 2,)) -# plus2x[1:ny - 1] = 1 / (2 * dx ** 2) -# # Associated RHS -# rhs[1:ny - 1] = rhs[1:ny - 1] + bcNeuArray[0] / dx -# -# # Check -# min2x = np.zeros((ny,)) -# min2x = np.tile(min2x, (nx - 2,)) -# min2x[-ny:] = 1 / (2 * dx ** 2) -# # Associated RHS -# rhs[-ny + 1:-1] = rhs[-ny + 1:-1] - bcNeuArray[1] / dx -# -# min1y = np.ones((ny,)) / (dy ** 2) -# min1y[0] = 0 # Dirichlet -# min1y[-2] = -2 / (dy ** 2) # One-side difference -# min1y[-1] = 0 # Not part of current column -# min1y = np.tile(min1y, (nx,))[0:-1] -# # Associated RHS -# rhs[1::ny] = rhs[1::ny] - bcDirArray[0] / (dy ** 2) -# -# plus1y = np.ones((ny,)) / (dy ** 2) -# plus1y[0] = -2 / (dy ** 2) # One-side difference -# plus1y[-2] = 0 # Dirichlet -# plus1y[-1] = 0 # Not part of current column -# plus1y = np.tile(plus1y, (nx,))[0:-1] -# # Associated RHS -# rhs[ny - 2::ny] = rhs[ny - 2::ny] - bcDirArray[1] / (dy ** 2) -# -# min2y = np.zeros((ny,)) -# min2y[1] = 1 / (dy ** 2) -# min2y = np.tile(min2y, (nx,))[0:-2] -# -# plus2y = np.zeros((ny,)) -# plus2y[0] = 1 / (dy ** 2) -# plus2y = np.tile(plus2y, (nx,))[0:-2] -# -# A = np.diag(mid) + np.diag(min1y, -1) + np.diag(min2y, -2) + np.diag(plus1y, 1) + np.diag(plus2y, 2) + np.diag( -# plus1x, ny) + np.diag( -# min1x, -ny) + \ -# np.diag(plus2x, ny * 2) + np.diag(min2x, -ny * 2) -# -# return A, rhs - - # Streamfunction operators def dxOpStream(): """Create discrete operator on 2d grid for dPsi / dx using central and one-side difference formulas. The one-sided @@ -396,4 +316,4 @@ def dlOpStreamMod(): mid[0:ny] = 1 mid[-ny:-1] = 1 - return sparse.csr_matrix(sparse.diags([mid, plus1y, min1y, plus1x, min1x], [0, 1, -1, ny, -ny])) + return sparse.csc_matrix(sparse.diags([mid, plus1y, min1y, plus1x, min1x], [0, 1, -1, ny, -ny]))