Simulation with Half Frap Initial conditions, prescribe Phi_u

5b7ca3a1 · MPI PKS Laptop · 4559fc85 · 5b7ca3a1
Commit 5b7ca3a1 authored 5 years ago by MPI PKS Laptop
--- a/Half_Frap_Phi_u.ipynb
+++ b/Half_Frap_Phi_u.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "#### READ ME ###\n",
+    "\n",
+    "################ When load this in Paraview !!!!!!!!!!!! => Work the first time then bug !!!!!!!!!!!\n",
+    "\n",
+    "from __future__ import print_function\n",
+    "from fenics import *\n",
+    "from mshr import *\n",
+    "from dolfin import *\n",
+    "import matplotlib.pyplot as plt\n",
+    "import time\n",
+    "from math import *\n",
+    "\n",
+    "set_log_level(40)\n",
+    "\n",
+    "\n",
+    "dt = 0.01\n",
+    "Radius_Droplet = 5\n",
+    "\n",
+    "\n",
+    "mesh = generate_mesh(Circle(Point(0,0),20.0),100)\n",
+    "\n",
+    "V = FunctionSpace(mesh, 'CG', 1)\n",
+    "Phi_u = Function(V)\n",
+    "v = TestFunction(V)\n",
+    "\n",
+    "\n",
+    "Phi_u0 = Function(V)\n",
+    "Phi_tot = Function(V)\n",
+    "\n",
+    "tol = 1E-14\n",
+    "\n",
+    "\n",
+    "Phi_u0_int = 0.01\n",
+    "Phi_u0_ext = 0.05\n",
+    "Phi_u0_notbleach = 0.35\n",
+    "\n",
+    "Phi_u0 = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet && atan2(x[1],x[0]) < 0  ? Phi_u0_int : sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet && atan2(x[1],x[0]) >= 0 ?  Phi_u0_notbleach : Phi_u0_ext' , degree=2,tol=tol,Phi_u0_int = Phi_u0_int,Phi_u0_ext=Phi_u0_ext, Phi_u0_notbleach = Phi_u0_notbleach, Radius_Droplet=Radius_Droplet),V)\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "Phi_tot_int = 0.35\n",
+    "Phi_tot_ext = 0.05\n",
+    "\n",
+    "Phi_tot = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet ? Phi_tot_int : Phi_tot_ext', degree=2,tol=tol,Phi_tot_int=Phi_tot_int,Phi_tot_ext=Phi_tot_ext,Radius_Droplet=Radius_Droplet),V)\n",
+    "\n",
+    "\n",
+    "u_D = Constant(0.05)\n",
+    "\n",
+    "def boundary(x, on_boundary):\n",
+    "    tol = 1E-14\n",
+    "    return on_boundary\n",
+    "\n",
+    "bc = DirichletBC(V, u_D, boundary)\n",
+    "\n",
+    "Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx\n",
+    "\n",
+    "\n",
+    "\n",
+    "vtkfile = File('Disk/Phi0.pvd')\n",
+    "vtkfile << Phi_u0\n",
+    "\n",
+    "\n",
+    "\n",
+    "t = 0\n",
+    "cFile = XDMFFile('test.xmdf')\n",
+    "cFile.write(Phi_u0, t)\n",
+    "\n",
+    "ti = time.time()\n",
+    "for i in range(200):\n",
+    "    solve(Form == 0, Phi_u, bc)\n",
+    "    assign(Phi_u0, Phi_u)\n",
+    "    t += dt\n",
+    "    cFile.write(Phi_u0, t)\n",
+    "    print(time.time() - ti)\n",
+    "\n",
+    "cFile.close()\n",
+    "\n",
+    "    \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "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.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
+%% Cell type:code id: tags:
+
+``` python
+
+#### READ ME ###
+
+################ When load this in Paraview !!!!!!!!!!!! => Work the first time then bug !!!!!!!!!!!
+
+from __future__ import print_function
+from fenics import *
+from mshr import *
+from dolfin import *
+import matplotlib.pyplot as plt
+import time
+from math import *
+
+set_log_level(40)
+
+
+dt = 0.01
+Radius_Droplet = 5
+
+
+mesh = generate_mesh(Circle(Point(0,0),20.0),100)
+
+V = FunctionSpace(mesh, 'CG', 1)
+Phi_u = Function(V)
+v = TestFunction(V)
+
+
+Phi_u0 = Function(V)
+Phi_tot = Function(V)
+
+tol = 1E-14
+
+
+Phi_u0_int = 0.01
+Phi_u0_ext = 0.05
+Phi_u0_notbleach = 0.35
+
+Phi_u0 = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet && atan2(x[1],x[0]) < 0  ? Phi_u0_int : sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet && atan2(x[1],x[0]) >= 0 ?  Phi_u0_notbleach : Phi_u0_ext' , degree=2,tol=tol,Phi_u0_int = Phi_u0_int,Phi_u0_ext=Phi_u0_ext, Phi_u0_notbleach = Phi_u0_notbleach, Radius_Droplet=Radius_Droplet),V)
+
+
+
+
+
+
+
+Phi_tot_int = 0.35
+Phi_tot_ext = 0.05
+
+Phi_tot = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet ? Phi_tot_int : Phi_tot_ext', degree=2,tol=tol,Phi_tot_int=Phi_tot_int,Phi_tot_ext=Phi_tot_ext,Radius_Droplet=Radius_Droplet),V)
+
+
+u_D = Constant(0.05)
+
+def boundary(x, on_boundary):
+    tol = 1E-14
+    return on_boundary
+
+bc = DirichletBC(V, u_D, boundary)
+
+Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx
+
+
+
+vtkfile = File('Disk/Phi0.pvd')
+vtkfile << Phi_u0
+
+
+
+t = 0
+cFile = XDMFFile('test.xmdf')
+cFile.write(Phi_u0, t)
+
+ti = time.time()
+for i in range(200):
+    solve(Form == 0, Phi_u, bc)
+    assign(Phi_u0, Phi_u)
+    t += dt
+    cFile.write(Phi_u0, t)
+    print(time.time() - ti)
+
+cFile.close()
+
+
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```