From be689856007ba31d92733693b5638efd3a807246 Mon Sep 17 00:00:00 2001 From: philippun Date: Thu, 22 Jun 2023 11:57:46 +0200 Subject: [PATCH] FTCS prototyping --- proto/FTCS.ipynb | 249 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 249 insertions(+) diff --git a/proto/FTCS.ipynb b/proto/FTCS.ipynb index e69de29..857b8e8 100644 --- a/proto/FTCS.ipynb +++ b/proto/FTCS.ipynb @@ -0,0 +1,249 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# imports\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "# environment\n", + "dim = 2\n", + "grid_size = {'x': 20, 'y': 20} # to start out\n", + "domain_size = {'x': 20, 'y': 20} # to start out\n", + "\n", + "delta_x = domain_size['x'] / grid_size['x']\n", + "delta_y = domain_size['y'] / grid_size['y'] \n", + "\n", + "alpha_x = np.ones((grid_size['x'],grid_size['y']))\n", + "alpha_y = np.ones((grid_size['x'],grid_size['y']))\n", + "\n", + "max_stable_time_step_x = delta_x**2 / (2 * np.max(alpha_x))\n", + "max_stable_time_step_y = delta_y**2 / (2 * np.max(alpha_y))\n", + "max_stable_time_step = max_stable_time_step_x if max_stable_time_step_x <= max_stable_time_step_y else max_stable_time_step_y\n", + "\n", + "C_t = np.zeros((grid_size['x'],grid_size['y']))" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "# initialization\n", + "C_t[1,1] = 1" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "def alpha_interblock(alpha1, alpha2, harmonic=False):\n", + " if not harmonic:\n", + " return 0.5 * (alpha1 + alpha2)\n", + " else:\n", + " return 2 / ((1/alpha1) + (1/alpha2))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "# simulating one time step\n", + "def simulate(C_t): \n", + " C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n", + "\n", + " # inner cells\n", + " for i in range(1, grid_size['x']-1):\n", + " for j in range(1, grid_size['y']-1):\n", + " C_t1[i,j] = C_t[i,j] \\\n", + " + max_stable_time_step/delta_x**2 * (alpha_interblock(alpha_x[i+1,j], alpha_x[i,j]) * C_t[i+1,j]\n", + " - (alpha_interblock(alpha_x[i+1,j], alpha_x[i,j]) + alpha_interblock(alpha_x[i-1,j], alpha_x[i,j])) * C_t[i,j] \n", + " + alpha_interblock(alpha_x[i-1,j], alpha_x[i,j]) * C_t[i-1,j]) \\\n", + " + max_stable_time_step/delta_y**2 * (alpha_interblock(alpha_y[i,j+1], alpha_y[i,j]) * C_t[i,j+1]\n", + " - (alpha_interblock(alpha_y[i,j+1], alpha_y[i,j]) + alpha_interblock(alpha_y[i,j-1], alpha_y[i,j])) * C_t[i,j] \n", + " + alpha_interblock(alpha_y[i,j-1], alpha_y[i,j]) * C_t[i,j-1])\n", + "\n", + " return C_t1" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 4.93210938e+02 -6.85390625e+02 5.61357422e+02\n", + " -3.21421875e+02 1.35512695e+02 -4.28554688e+01 1.01855469e+01\n", + " -1.79687500e+00 2.28515625e-01 -1.95312500e-02 9.76562500e-04\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 -6.85390625e+02 9.41660156e+02 -7.56210938e+02\n", + " 4.20468750e+02 -1.70039062e+02 5.07128906e+01 -1.10742188e+01\n", + " 1.71875000e+00 -1.75781250e-01 9.76562500e-03 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 5.61357422e+02 -7.56210938e+02 5.86933594e+02\n", + " -3.10078125e+02 1.16542969e+02 -3.12890625e+01 5.84472656e+00\n", + " -7.03125000e-01 4.39453125e-02 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 -3.21421875e+02 4.20468750e+02 -3.10078125e+02\n", + " 1.51593750e+02 -5.08593750e+01 1.15312500e+01 -1.64062500e+00\n", + " 1.17187500e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 1.35512695e+02 -1.70039062e+02 1.16542969e+02\n", + " -5.08593750e+01 1.43554688e+01 -2.46093750e+00 2.05078125e-01\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 -4.28554688e+01 5.07128906e+01 -3.12890625e+01\n", + " 1.15312500e+01 -2.46093750e+00 2.46093750e-01 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 1.01855469e+01 -1.10742188e+01 5.84472656e+00\n", + " -1.64062500e+00 2.05078125e-01 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 -1.79687500e+00 1.71875000e+00 -7.03125000e-01\n", + " 1.17187500e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 2.28515625e-01 -1.75781250e-01 4.39453125e-02\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 -1.95312500e-02 9.76562500e-03 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 9.76562500e-04 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", + " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", + " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n" + ] + } + ], + "source": [ + "for i in range(10):\n", + " C_t = simulate(C_t)\n", + "print(C_t)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.5\n" + ] + } + ], + "source": [ + "print(max_stable_time_step)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.1" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +}