From 37d16d13a495d37d3577cafae41f06d417a0f776 Mon Sep 17 00:00:00 2001 From: Lars Gebraad Date: Thu, 22 Mar 2018 11:41:38 +0100 Subject: [PATCH] Removed legacy code. --- alt_main.py | 80 +++++++++--- boundary.py | 89 ------------- kernels.py | 288 ------------------------------------------- main.py | 93 -------------- nonlinsys.py | 15 --- onedimensionalode.py | 43 ------- 6 files changed, 60 insertions(+), 548 deletions(-) delete mode 100644 boundary.py delete mode 100644 kernels.py delete mode 100755 main.py delete mode 100644 nonlinsys.py delete mode 100644 onedimensionalode.py diff --git a/alt_main.py b/alt_main.py index 74a54fa..64b868f 100644 --- a/alt_main.py +++ b/alt_main.py @@ -2,36 +2,38 @@ import matplotlib.pyplot as plt import numpy as np import scipy.sparse as sparse -import scipy.linalg as sparselinalg 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 -dt = 0.05 +# --- 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) -# --- Setup --- # -config.ny = 80 -config.nx = config.ny * 2 -config.dy = 1 -config.dx = 1 -# 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 -# start with linear profile for temperature 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) @@ -41,8 +43,8 @@ dxOpPsi = op.dxOpStream() dlOpPsi = op.dlOpStreamMod() -dlOpTemp, rhsDlOpTemp = op.dlOpTemp(bcDirArray=tempDirichlet, bcNeuArray=[0, 0]) -dxOpTemp, rhsDxOpTemp = op.dxOpTemp(bcNeuArray=[0, 0]) +dlOpTemp, rhsDlOpTemp = op.dlOpTemp(bcDirArray=tempDirichlet, bcNeuArray=tempNeumann) +dxOpTemp, rhsDxOpTemp = op.dxOpTemp(bcNeuArray=tempNeumann) dyOpTemp, rhsDyOpTemp = op.dyOpTemp(bcDirArray=tempDirichlet) Tnplus1 = startT @@ -52,22 +54,60 @@ toKeep[0:ny] = 0 toKeep[-ny:] = 0 toKeep[ny::ny] = 0 -toKeep[ny-1::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(100): +for it in range(500): Tn = Tnplus1 Psin = Psinplus1 - C, rhsC = matrices.constructC(dxOpTemp=dxOpTemp, rhsDxOpTemp=rhsDxOpTemp, dyOpTemp=dyOpTemp, rhsDyOpTemp=rhsDyOpTemp, + 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) - Tnplus1.shape = (nx*ny,1) - Psinplus1 = sparse.linalg.spsolve(dlOpPsi,-sqrtRa * toKeep @ (dxOpTemp @ Tnplus1 - rhsDxOpTemp)) - Psinplus1.shape = (nx*ny,1) + dt * rhsC + (sparse.eye(nx * ny) - (dt / 2) * C) @ Tn, use_umfpack=True) + + # Enforcing column vector + Tnplus1.shape = (nx * ny, 1) -plt.imshow(grid.vecToField(Tnplus1)) + # 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/boundary.py b/boundary.py deleted file mode 100644 index cb6fe1e..0000000 --- a/boundary.py +++ /dev/null @@ -1,89 +0,0 @@ -import numpy as np -import scipy.sparse as sparse - -import config - - -class BoundaryValue(): - def load(self): - self.ny = config.ny - self.nx = config.nx - self.dy = config.dy - self.dx = config.dx - - def __init__(self, x0=None, xend=None, y0=None, yend=None): - """ - Setting Dirichlet Boundary conditions. The hierarchy of these is x0, xend, y0, yend - This means that if both boundaries are not equal, the first in this list gets priority. - :param x0: value at x0, default is None - :param xend: value at xend, default is None - :param y0: value at y0, default is None - :param yend: value at yend, default is None - """ - self.load() - - row = np.array([]) - col = np.array([]) - data = np.array([]) - self.fixedVal = np.array([0, 0, 0, 0]) - - if (x0 != None): - row = np.concatenate((row, np.arange(self.ny))) - col = np.concatenate((col, np.zeros(self.ny))) - data = np.concatenate((data, np.ones(self.ny))) - self.fixedVal[0] = x0 - if (xend != None): - row = np.concatenate((row, np.arange(self.ny))) - col = np.concatenate((col, (self.nx - 1) * np.ones(self.ny))) - data = np.concatenate((data, 2 * np.ones(self.ny))) - self.fixedVal[1] = xend - if (y0 != None): - y0start = 0 - numy0 = self.nx - # x0 boundary override - if (x0 != None): - numy0 = numy0 - 1 - y0start = y0start + 1 - # xend boundary override - if (xend != None): - numy0 = numy0 - 1 - - row = np.concatenate((row, np.zeros(numy0))) - col = np.concatenate((col, y0start + np.arange(numy0))) - data = np.concatenate((data, 3 * np.ones(numy0))) - self.fixedVal[2] = y0 - if (yend != None): - yendstart = 0 - numyend = self.nx - # x0 boundary override - if (x0 != None): - numyend = numyend - 1 - yendstart = yendstart + 1 - # xend boundary override - if (xend != None): - numyend = numyend - 1 - - row = np.concatenate((row, (self.ny - 1) * np.ones(numyend))) - col = np.concatenate((col, yendstart + np.arange(numyend))) - data = np.concatenate((data, 4 * np.ones(numyend))) - self.fixedVal[3] = yend - self.SparseIndices = sparse.dok_matrix(sparse.coo_matrix((data, (row, col)), shape=(self.ny, self.nx)), - dtype='int') - self.Values = [x0, xend, y0, yend] - # print('Sparsity: ', self.DirichletSparseMat.getnnz() / (self.nx * self.ny)) - - -class Dirichlet(BoundaryValue): - def __init__(self, x0=None, xend=None, y0=None, yend=None): - super().__init__(x0, xend, y0, yend) - self.type = "Dirichlet" - - -class Neumann(BoundaryValue): - def __init__(self, x0=None, xend=None, y0=None, yend=None): - if (y0 != None or yend != None): - # No need in this project - raise (NotImplementedError) - - super().__init__(x0, xend, y0, yend) - self.type = "Neumann" diff --git a/kernels.py b/kernels.py deleted file mode 100644 index ff0e93e..0000000 --- a/kernels.py +++ /dev/null @@ -1,288 +0,0 @@ -import scipy.sparse as sparse - -import config -import discretizationUtils as grid - - -def makeSparseLaplaceKernel(dirichlet=None, neumann=None): - ny = config.ny - nx = config.nx - dy = config.dy - dx = config.dx - Laplace = sparse.dok_matrix((nx * ny, nx * ny), dtype='float') - rhsLaplace = sparse.dok_matrix((nx * ny, 1), dtype='float') - - for index, row in enumerate(Laplace): - if (dirichlet.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - # print('Dirichlet value at pos %i %i' % (grid.index1to2y(index), grid.index1to2x(index))) - Laplace[index, index] = 1 - rhsLaplace[index] = dirichlet.fixedVal[ - dirichlet.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] - 1] - elif (neumann != None and neumann.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - # print('Neumann value at pos %i %i' % (grid.index1to2y(index), grid.index1to2x(index))) - if (grid.index1to2x(index) == 0 or grid.index1to2x(index) == nx - 1): - # Neumann in x-direction - Laplace[index, index] = -(2.0 / dx ** 2 + 2.0 / dy ** 2) - - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - # This implements a Laplace operator on the defined grid - try: - Laplace[index, grid.index2to1(curY - 1, curX)] = 1.0 / dy ** 2 - except (IndexError): - pass - - try: - Laplace[index, grid.index2to1(curY + 1, curX)] = 1.0 / dy ** 2 - except (IndexError): - pass - - # Which point in the x-sense don't we know? - if (grid.index1to2x(index) == 0): - try: - Laplace[index, grid.index2to1(curY, curX + 1)] = 2 / dx ** 2 - rhsLaplace[index] = - 2 * dx * neumann.fixedVal[ - neumann.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] - 1] - except (IndexError): - pass - else: - try: - Laplace[index, grid.index2to1(curY, curX - 1)] = 2 / dx ** 2 - rhsLaplace[index] = 2 * dx * neumann.fixedVal[ - neumann.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] - 1] - except (IndexError): - pass - else: - # Neumann in y-direction - raise (NotImplementedError) - else: - # print('No Dirichlet value at pos %i %i' % (grid.index1to2y(index), grid.index1to2x(index))) - Laplace[index, index] = -(2.0 / dx ** 2 + 2.0 / dy ** 2) - - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - # This implements a Laplace operator on the defined grid - try: - Laplace[index, grid.index2to1(curY - 1, curX)] = 1.0 / dy ** 2 - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - try: - Laplace[index, grid.index2to1(curY + 1, curX)] = 1.0 / dy ** 2 - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - try: - Laplace[index, grid.index2to1(curY, curX - 1)] = 1.0 / dx ** 2 - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - try: - Laplace[index, grid.index2to1(curY, curX + 1)] = 1.0 / dx ** 2 - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - return Laplace, rhsLaplace - - -def makeSparsedxKernel(dirichlet=None, neumann=None): - ny = config.ny - nx = config.nx - dy = config.dy - dx = config.dx - dxMat = sparse.dok_matrix((nx * ny, nx * ny), dtype='float') - rhsDx = sparse.dok_matrix((nx * ny, 1), dtype='float') - - for index, row in enumerate(dxMat): - if (neumann != None and neumann.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - # print('Neumann value at pos %i %i' % (grid.index1to2y(index), grid.index1to2x(index))) - if (grid.index1to2x(index) == 0 or grid.index1to2x(index) == nx - 1): - # Neumann in x-direction - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - rhsDx[index] = - neumann.fixedVal[neumann.SparseIndices[curY, curX] - 1] - - else: - # Neumann in y-direction - raise (NotImplementedError) - elif (grid.index1to2x(index) == 0): - # print('Neumann condition absent, needed for central difference at boundary. Falling back to one sided Euler...') - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - # Left point, let's try to implement a dirichlet value instead using one-sided euler - if (dirichlet != None and dirichlet.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - try: - rhsDx[index] = dirichlet.fixedVal[dirichlet.SparseIndices[curY, curX] - 1] / (dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - else: - try: - dxMat[index, grid.index2to1(curY, curX)] = -1.0 / (dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - # Right point - try: - dxMat[index, grid.index2to1(curY, curX + 1)] = 1.0 / (dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - elif (grid.index1to2x(index) == nx - 1): - # print('Neumann condition absent, needed for central difference at boundary. Falling back to one sided Euler...') - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - try: - dxMat[index, grid.index2to1(curY, curX - 1)] = -1.0 / (dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - if (dirichlet != None and dirichlet.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - try: - rhsDx[index] = - dirichlet.fixedVal[dirichlet.SparseIndices[curY, curX] - 1] / (dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - else: - try: - dxMat[index, grid.index2to1(curY, curX)] = 1.0 / (dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - else: - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - # Left point - try: - dxMat[index, grid.index2to1(curY, curX - 1)] = -1.0 / (2.0 * dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - # Right point - try: - dxMat[index, grid.index2to1(curY, curX + 1)] = 1.0 / (2.0 * dx) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - return dxMat, rhsDx - - -def makeSparsedyKernel(dirichlet=None, neumann=None): - ny = config.ny - nx = config.nx - dy = config.dy - dx = config.dx - dyMat = sparse.dok_matrix((nx * ny, nx * ny), dtype='float') - rhsDy = sparse.dok_matrix((nx * ny, 1), dtype='float') - - for index, row in enumerate(dyMat): - if (neumann != None and neumann.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - # print('Neumann value at pos %i %i' % (grid.index1to2y(index), grid.index1to2x(index))) - if (grid.index1to2y(index) == 0 or grid.index1to2y(index) == ny - 1): - # Neumann in y-direction - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - rhsDy[index] = - neumann.fixedVal[neumann.SparseIndices[curY, curX] - 1] - - # else: - # Neumann in x-direction - # raise (NotImplementedError) - - elif (grid.index1to2y(index) == 0): - # print( - # 'Neumann condition absent, needed for central difference at boundary. Falling back to one sided Euler...') - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - # Left point, let's try to implement a dirichlet value instead using one-sided euler - if (dirichlet != None and dirichlet.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - try: - rhsDy[index] = dirichlet.fixedVal[dirichlet.SparseIndices[curY, curX] - 1] / (dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - else: - try: - dyMat[index, grid.index2to1(curY, curX)] = -1.0 / (dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - # Right point - try: - dyMat[index, grid.index2to1(curY + 1, curX)] = 1.0 / (dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - elif (grid.index1to2y(index) == ny - 1): - # print( - # 'Neumann condition absent, needed for central difference at boundary. Falling back to one sided Euler...') - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - try: - dyMat[index, grid.index2to1(curY - 1, curX)] = -1.0 / (dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - if (dirichlet != None and dirichlet.SparseIndices[grid.index1to2y(index), grid.index1to2x(index)] > 0): - try: - rhsDy[index] = - dirichlet.fixedVal[dirichlet.SparseIndices[curY, curX] - 1] / (dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - else: - try: - dyMat[index, grid.index2to1(curY, curX)] = 1.0 / (dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - else: - curX = grid.index1to2x(index) - curY = grid.index1to2y(index) - - # Left point - try: - dyMat[index, grid.index2to1(curY - 1, curX)] = -1.0 / (2.0 * dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - # Right point - try: - dyMat[index, grid.index2to1(curY + 1, curX)] = 1.0 / (2.0 * dy) - except (IndexError): - # TODO make relevant exception, however this should never occur... - pass - - return dyMat, rhsDy - - -def plotKernel(mat, vec): - import matplotlib.pyplot as plt - f, axs = plt.subplots(1, 2) - f.set_figheight(10) - f.set_figwidth(10) - im1 = axs[0].matshow(mat.todense().A) - # plt.colorbar(im1,cax=axs[0]) - im2 = axs[1].matshow(vec.todense().A, aspect=0.1) - # plt.colorbar(im2,cax=axs[1]) - plt.show() diff --git a/main.py b/main.py deleted file mode 100755 index 69333ce..0000000 --- a/main.py +++ /dev/null @@ -1,93 +0,0 @@ -#!/usr/bin/python3.6 -import matplotlib.pyplot as plt -import numpy as np - -import boundary -import config -import discretizationUtils as grid -import kernels - -config.ny = 11 -config.nx = config.ny * 2 -config.dy = 1 -config.dx = 1 - -ny = config.ny -nx = config.nx -dy = config.dy -dx = config.dx - -sqrt_ra = 7 - -Tbottom = 1 -Ttop = 0 - -# Initial conditions -startT = np.expand_dims(np.linspace(0, 1, ny), 1).repeat(nx, axis=1) -startT = np.expand_dims(grid.fieldToVec(startT), 1) - -startPsi = np.expand_dims(np.zeros(ny), 1).repeat(nx, axis=1) -startPsi = np.expand_dims(grid.fieldToVec(startPsi), 1) - -## Definition of boundary value problem -dirTemp = boundary.Dirichlet(y0=Tbottom, yend=Ttop) -neuTemp = boundary.Neumann(xend=0, x0=0) -dirPsi = boundary.Dirichlet(x0=0, y0=0, xend=0, yend=0) - -# Define matrices needed for every derivative in systems of algebraic equations using programmed kernels and bc's -laplaceMatPsi, laplaceVecPsi = kernels.makeSparseLaplaceKernel(dirichlet=dirPsi, neumann=None) -laplaceMatT, laplaceVecT = kernels.makeSparseLaplaceKernel(dirichlet=dirTemp, neumann=neuTemp) - -dxMatPsi, dxVecPsi = kernels.makeSparsedxKernel(dirichlet=dirPsi, neumann=None) -dxMatT, dxVecT = kernels.makeSparsedxKernel(dirichlet=dirTemp, neumann=neuTemp) - -dyMatPsi, dyVecPsi = kernels.makeSparsedyKernel(dirichlet=dirPsi, neumann=None) -dyMatT, dyVecT = kernels.makeSparsedyKernel(dirichlet=dirTemp, neumann=neuTemp) - - -# Define residual -def toMinimize(input): - global sqrt_ra - global solution_n_min_1 - - psi_p, t_p = np.split(solution_n_min_1, 2) - psi, t = np.split(input, 2) - - a = laplaceMatPsi @ psi + sqrt_ra * dxMatT @ t - (laplaceVecPsi + sqrt_ra * dxVecT) - b = sqrt_ra * ( - np.multiply(dyMatPsi @ psi, dxMatT @ t) - np.multiply(dxMatPsi @ psi, dyMatT @ t) - ) - laplaceMatPsi @ psi - (sqrt_ra * (dxVecT.multiply(dyVecPsi) - dyVecT.multiply(dxVecPsi)) + laplaceVecPsi) - - # a_prev = laplaceMatPsi @ psi_p + sqrt_ra * dxMatT @ t_p - (laplaceVecPsi + sqrt_ra * dxVecT) - - b_prev = sqrt_ra * ( - np.multiply(dyMatPsi @ psi_p, dxMatT @ t_p) - np.multiply(dxMatPsi @ psi_p, dyMatT @ t_p) - ) - laplaceMatPsi @ psi_p - (sqrt_ra * (dxVecT.multiply(dyVecPsi) - dyVecT.multiply(dxVecPsi)) + laplaceVecPsi) - - return np.concatenate((a, b - b_prev)) - - -currentPsi = startPsi -currentT = startT - -currentSol = np.concatenate((currentPsi, currentT)) - - -dt = 0.5 -dT = laplaceMatT @ currentT + laplaceVecT - sqrt_ra * ( - np.multiply( (dyMatPsi @ currentPsi + dyVecPsi), (dxMatT @ currentT + dxVecT) ) - - - np.multiply( (dxMatPsi @ currentPsi + dxVecPsi), (dyMatT @ currentT + dyVecT) ) -) - -nextT = currentT + dT * dt - -currentT[2] = 2 - -plt.matshow(grid.vecToField(currentT)) -plt.colorbar() -plt.show() - -plt.matshow(grid.vecToField(dyMatT @ currentT + dyVecT[::-1])) -plt.colorbar() -plt.show() \ No newline at end of file diff --git a/nonlinsys.py b/nonlinsys.py deleted file mode 100644 index ad90f9c..0000000 --- a/nonlinsys.py +++ /dev/null @@ -1,15 +0,0 @@ -import numpy as np -from scipy.optimize import fsolve - - -def equations(p): - x = p - np.arange(0, 10) - return x - - -x0 = np.arange(0, 10) + 3 -y0 = np.arange(0, 10) + 3 - -xsol = fsolve(equations, x0) - -print(equations((xsol))) diff --git a/onedimensionalode.py b/onedimensionalode.py deleted file mode 100644 index 2aac90c..0000000 --- a/onedimensionalode.py +++ /dev/null @@ -1,43 +0,0 @@ -#!/usr/bin/python3.6 -import matplotlib.pyplot as plt -import numpy as np -import scipy.optimize as optimize - -import boundary -import config -import discretizationUtils as grid -import kernels - -config.ny = 10 -config.nx = 10 -config.dy = 0.1 -config.dx = 0.1 - -ny = config.ny -nx = config.nx -dy = config.dy -dx = config.dx - -## Definition of boundary value problem -neuPsi = None -dirPsi = boundary.Dirichlet(y0=0) - -# Define matrices needed for every derivative in systems of algebraic equations using programmed kernels and bc's -dyMatPsi, dyVecPsi = kernels.makeSparsedyKernel(dirichlet=dirPsi, neumann=neuPsi) - -kernels.plotKernel(dyMatPsi, dyVecPsi) - - -# Define residual -def toMinimize(solution): - return (dyMatPsi @ solution - dyVecPsi - 1) - - -startT = np.expand_dims(np.zeros(ny), 1).repeat(nx, axis=1) - -sol = optimize.broyden2(toMinimize, np.vstack( - (np.expand_dims(grid.fieldToVec(startT), 1), np.expand_dims(grid.fieldToVec(startT), 1))), f_tol=6e-6) - -plt.matshow(grid.vecToField(sol)) -plt.colorbar() -plt.show()