diff --git a/proto/BTCS.ipynb b/proto/BTCS.ipynb index 7659196..dc8ff0b 100644 --- a/proto/BTCS.ipynb +++ b/proto/BTCS.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 4, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -68,7 +68,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -82,7 +82,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -92,88 +92,145 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 38, "metadata": {}, "outputs": [], "source": [ + "# CONSTANT\n", + "\n", "# creates the coeffiecient matrix for a given row\n", "def create_coeff_matrix(alpha_x, row_index, s_x, transpose=False):\n", + " \n", + "\n", + " r = row_index\n", " if (transpose):\n", " alpha_x = alpha_x.transpose()\n", "\n", - " cm = np.zeros((alpha_x.shape[0]+2,alpha_x.shape[1]+2))\n", + " # this is always a square matrix -> column by column, so \n", + " cm = np.zeros((alpha_x.shape[1],alpha_x.shape[1]))\n", "\n", - " offset = 2\n", - " r = row_index\n", - " for i in range(2, cm.shape[0]-2):\n", - " j = offset\n", - " cm[i,j-1] = -s_x * alpha_interblock(alpha_x[r, j-1-1], alpha_x[r, j-1])\n", - " cm[i,j] = 1 + s_x * (alpha_interblock(alpha_x[r,j-1], alpha_x[r, j+1-1]) + alpha_interblock(alpha_x[r, j-1], alpha_x[r, j-1-1]))\n", - " cm[i,j+1] = -s_x * alpha_interblock(alpha_x[r, j-1], alpha_x[r, j+1-1])\n", + " # top left\n", + " cm[0,0] = 1 + s_x * (alpha_interblock(alpha_x[r,0], alpha_x[r,1]) + 2 * alpha_x[r,0]) # alpha_x[0,0]\n", + " cm[0,1] = -s_x * alpha_interblock(alpha_x[r, 0], alpha_x[r, 1])\n", "\n", - " offset += 1\n", + " # inner cells\n", + " for i_cm in range(1, cm.shape[0]-1):\n", + " j_alpha = i_cm - 1\n", + " cm[i_cm,i_cm-1] = -s_x * alpha_interblock(alpha_x[r, j_alpha-1], alpha_x[r, j_alpha])\n", + " cm[i_cm,i_cm] = 1 + s_x * (alpha_interblock(alpha_x[r,j_alpha], alpha_x[r, j_alpha+1]) + alpha_interblock(alpha_x[r, j_alpha], alpha_x[r, j_alpha-1]))\n", + " cm[i_cm,i_cm+1] = -s_x * alpha_interblock(alpha_x[r, j_alpha], alpha_x[r, j_alpha+1])\n", + "\n", + "\n", + " # bottom right\n", + " n = cm.shape[0]-1\n", + " cm[n,n-1] = s_x * alpha_interblock(alpha_x[r, n-1], alpha_x[r, n])\n", + " cm[n,n] = 1 + s_x * (-alpha_interblock(alpha_x[r, n], alpha_x[r, n-1]) + 2 * alpha_x[n,n]) # alpha_x[n,n]?\n", "\n", " return cm\n", "\n", "# creates the solution vector for a given row\n", - "def create_solution_vector(concentrations, alpha_y, row_index, s_x, transpose=False):\n", - " pass\n", + "# s_x, alpha_x, row_index?\n", + "def create_solution_vector(concentrations, alpha_y, alpha_x, col_index, row_index, s_y, s_x, bc_left, bc_right, transpose=False):\n", + "\n", + " # \n", + " sv = np.zeros(alpha_y.shape[1])\n", + "\n", + " # inner\n", + " for i_sv in range(0, sv.shape[0]):\n", + " i_alpha = i_sv-1\n", + " sv[i_sv] = s_y * alpha_interblock(alpha_y[i_alpha+1,col_index],alpha_y[i_alpha,col_index])*concentrations[i_alpha+1,col_index] \\\n", + " + (1 - s_y *(alpha_interblock(alpha_y[i_alpha+1,col_index], alpha_y[i_alpha,col_index])) \\\n", + " + alpha_interblock(alpha_y[i_alpha,col_index], alpha_y[i_alpha-1,col_index]))*concentrations[i_alpha,col_index] \\\n", + " + s_y * alpha_interblock(alpha_y[i_alpha,col_index], alpha_y[i_alpha-1,col_index])*concentrations[i_alpha-1,col_index]\n", + "\n", + " # first\n", + " sv[0] += 2 * s_x * alpha_x[row_index,0] * bc_left # alpha_x[row_index,0]?\n", + "\n", + " # last\n", + " n = sv.shape[0]-1\n", + " sv[n] += 2 * s_x * alpha_x[row_index,n] * bc_right\n", + "\n", + " return sv\n", "\n", "\n", "def solve_row(A, b):\n", - " # use existing implementation\n", - " pass" + " return np.linalg.solve(A,b)\n", + "\n", + "def calc_sweep(concentrations, alpha_x, alpha_y, s_y, s_x, bc_left, bc_right, transpose=False):\n", + " for row in range(concentrations.shape[0]):\n", + " A = create_coeff_matrix(alpha_x, row, s_x, transpose)\n", + " b = create_solution_vector(concentrations, alpha_y, alpha_x, )\n", + "\n", + "def calc_round():\n", + " calc_sweep()\n", + " calc_sweep()\n" ] }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "[[ 0. 0. 0. 0. 0. 0.\n", - " 0. ]\n", - " [ 0. 0. 0. 0. 0. 0.\n", - " 0. ]\n", - " [ 0. -1.28906303 3.7948142 -1.50575117 0. 0.\n", - " 0. ]\n", - " [ 0. 0. -1.50575117 4.6735985 -2.16784733 0.\n", - " 0. ]\n", - " [ 0. 0. 0. -2.16784733 5.48692172 -2.31907439\n", - " 0. ]\n", - " [ 0. 0. 0. 0. 0. 0.\n", - " 0. ]\n", - " [ 0. 0. 0. 0. 0. 0.\n", - " 0. ]]\n" + "Test A:\n", + "[[ 2.70902018 -0.77802779 0. 0. 0. ]\n", + " [-0.34610519 2.12413297 -0.77802779 0. 0. ]\n", + " [ 0. -0.77802779 3.74324914 -1.96522135 0. ]\n", + " [ 0. 0. -1.96522135 4.21016215 -1.2449408 ]\n", + " [ 0. 0. 0. 0.43089455 5.13559523]]\n", + "\n", + "TestB:\n", + "[2.45811954 1.70379428 2.70189936 0.92506681 4.0009449 ]\n", + "\n", + "Result\n", + "[1.38151174 1.65087116 1.70267339 1.21472714 0.67714168]\n" ] } ], "source": [ - "result = create_coeff_matrix(np.random.random_sample((5,5)), 0, 3)\n", + "row = 5\n", + "col = 5\n", + "\n", + "testConcentrations = np.random.random((row,col))\n", + "testAlpha = np.random.random((row,col))\n", + "testA = create_coeff_matrix(testAlpha, 1, 3)\n", + "print(\"Test A:\")\n", + "print(testA)\n", + "# print(test.shape[0])\n", + "testb = create_solution_vector(testConcentrations,testAlpha,testAlpha,1,1,3,3,1,1)\n", + "print(\"\\nTestB:\")\n", + "print(testb)\n", + "\n", + "# needed to solve: TODO calc boundary in functions\n", + "# testA[0,0] = 1\n", + "# testA[4,4] = 1\n", + "# testb[0] = 1\n", + "# test[4] = 1\n", + "\n", + "result = solve_row(testA, testb)\n", + "print(\"\\nResult\")\n", "print(result)" ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "2\n", - "3\n", - "4\n" + "3\n" ] } ], "source": [ - "for i in range(2,5):\n", - " print(i)" + "test2 = np.random.random_sample(3)\n", + "print(test2.shape[0])" ] }, {