{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Level-set method\n",
"***"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"The level-set method is a method designed for modeling material interaction. The main idea is based on a function $l$, which satisfies the following.\n",
"\n",
"$$ l(x) \\begin{cases}\n",
" > 0 \\text{ if } x \\in \\Omega_1 \\\\\n",
" < 0 \\text{ in } x \\in \\Omega_2\n",
" \\end{cases}. $$\n",
"\n",
"We will focus only on the interaction between two incompressible Navier-Stokes equations with different\n",
"constant parameters $\\rho_i, \\mu_i \\in \\Omega_i$ for $i \\in \\{ 1, 2 \\}$. The function $l$ allows us to define\n",
"viscosity $\\mu$ and density $\\rho$ in the whole domain $\\Omega$ as a single\n",
"function\n",
"\\begin{equation*}\n",
" \\mu(l) = \\frac{1}{2}(\\text{sign}(l) + 1)\\mu_1 + \\frac{1}{2}(\\text{sign}(l) - 1) \\mu_2\n",
"\\end{equation*}\n",
"and\n",
"\\begin{equation*}\n",
" \\rho(l) = \\frac{1}{2}(\\text{sign}(l) + 1)\\rho_1 + \\frac{1}{2}(\\text{sign}(l) - 1) \\rho_2.\n",
"\\end{equation*}\n",
"Because we would like to solve the evolution problem, the $\\Omega_i$ and the function $l$\n",
"have to depend on time. We want $l$ to be constant along streamlines for the prescribed velocity field $v$.\n",
"This condition is satisfied if\n",
"\n",
"\\begin{equation*}\n",
"\\partial_t l(x, t) + \\text{div}(l(x, t) v(x, t)) = 0.\n",
"\\end{equation*}\n",
"\n",
"We can now formulate the whole system, firstly in the strong sense.\n",
"$$ \\rho(l) \\left( \\partial_t v + (v \\cdot \\nabla) v \\right) = \\text{div}\\left(\\mathbb{T}(\\mu(l), \\nabla v)\\right) + \\rho(l)g,$$\n",
"\n",
"$$ \\mathbb{T}(\\mu(l), \\nabla v) = \\mu(l)\\left( \\nabla v + (\\nabla v)^T \\right) - p\\mathbb{I}, $$\n",
"\n",
"$$ \\partial_t l + \\text{div}(l v) = 0, $$\n",
"\n",
"$$ \\text{div}(v) = 0. $$"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rayleigh-Taylor Example\n",
"Let us assume a rectangular domain with two fluids as it is desctibed in the figure below. We will cosidere Navier-Stokes fluid for both of them.\n",
"The fluids will stick on the boundary of the domain, so $v = 0$ at $\\partial \\Omega$.\n",
"\n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"mesh = df.RectangleMesh(\n",
" df.Point(0.0, -0.5), df.Point(0.25, 0.5), 20, 80, 'crossed'\n",
")\n",
"df.plot(mesh)\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We define the initial conditions for velocity and levelset. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"class InitialCondition(df.UserExpression):\n",
"\n",
" def __init__(self, **kwargs):\n",
" self.center = [0.125, 0.25]\n",
" self.radius = df.sqrt(0.125**2 + 0.25**2)\n",
" super().__init__(**kwargs)\n",
" \n",
" def r(self, x: list) -> float:\n",
" return df.sqrt((x[0] - self.center[0])**2 + (x[1] - self.center[1])**2)\n",
"\n",
" def eval(self, values, x):\n",
" values[0] = 0.0 # v_x\n",
" values[1] = 0.0 # v_y\n",
" values[2] = 0.0 # p\n",
" values[3] = self.r(x) - self.radius # l\n",
"\n",
" def value_shape(self):\n",
" return (4, )\n",
"# create instance of the class.\n",
"# initial_conditions = InitialCondition()\n",
"\n",
"\n",
"class InitialConditionSigmoid(df.UserExpression):\n",
"\n",
" def __init__(self, **kwargs):\n",
" self.center = [0.125, 0.25]\n",
" self.radius = df.sqrt(0.125**2 + 0.25**2)\n",
" self.eps = 0.001 \n",
" super().__init__(**kwargs)\n",
" \n",
" def r(self, x: list) -> float:\n",
" return df.sqrt((x[0] - self.center[0])**2 + (x[1] - self.center[1])**2)\n",
"\n",
" def eval(self, values, x):\n",
" values[0] = 0.0 # v_x\n",
" values[1] = 0.0 # v_y\n",
" values[2] = 0.0 # p\n",
" values[3] = 1 / (1 + df.exp(min((self.radius - self.r(x) ) / self.eps, 10)))\n",
" # values[3] = 1 / (df.exp( (self.r(x) - self.radius) / self.eps ))\n",
"\n",
" def value_shape(self):\n",
" return (4, )\n",
"\n",
"initial_conditions = InitialConditionSigmoid()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we define approximation of the signum function as\n",
"\\begin{equation}\n",
" \\text{sign}(l) = \\frac{l}{\\sqrt(l^2 + \\varepsilon^2\\nabla l \\cdot \\nabla l)} \n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# def sign(q: df.Function, eps: float):\n",
"# return q / df.sqrt(q * q + eps * eps * df.inner(df.grad(q), df.grad(q)))\n",
"\n",
"def sign(q: df.Function, eps: float):\n",
" return 2 * (q - 0.5) "
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Further we create the function spaces. We use \"CG\" for each subspace."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"elements = [\n",
" df.VectorElement(\"CG\", mesh.ufl_cell(), 2), # velocity\n",
" df.FiniteElement(\"CG\", mesh.ufl_cell(), 1), # pressure\n",
" df.FiniteElement(\"CG\", mesh.ufl_cell(), 1) # levelset\n",
"]\n",
"\n",
"function_space = df.FunctionSpace(\n",
" mesh, df.MixedElement(elements)\n",
")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we define all the constants. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"material_params = {\n",
" \"mu1\": 1.0,\n",
" \"mu2\": 1.0,\n",
" \"rho1\": 500,\n",
" \"rho2\": 1000,\n",
"}\n",
"eps = 1e-4\n",
"dt = 0.02\n",
"t_start = 0.0\n",
"t_end = 1.0\n",
"g = df.Constant((0.0, -10.0)) # gravity field"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Further, we create the boundary conditions for velocity. In particular we \n",
"require that $v = 0$ on the whole boundary of the domain."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"bcs = [\n",
" df.DirichletBC(\n",
" function_space.sub(0), # v\n",
" df.Constant((0.0, 0.0)),\n",
" \"on_boundary\"\n",
" )\n",
"]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we define the funcions on our space and interpolate the initial condition on the space."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"w = df.Function(function_space) # unknown\n",
"w0 = df.Function(function_space) # from previous step\n",
"phi = df.TestFunction(function_space) # test function\n",
"w0.assign(df.interpolate(initial_conditions, function_space))\n",
"w.assign(w0)\n",
"\n",
"# Split functions\n",
"v, p, l = df.split(w)\n",
"v0, p0, l0 = df.split(w0)\n",
"phi_v, phi_p, phi_l = df.split(phi)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We will need to formulate the functions $\\rho$ and $\\mu$. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def rho(params: dict, l: df.Function, eps: float, sign):\n",
" return (\n",
" params[\"rho1\"] * 0.5* (1.0 + sign(l, eps))\n",
" + params[\"rho2\"] * 0.5 * (1.0 - sign(l, eps))\n",
" )\n",
"\n",
"\n",
"def mu(params: dict, l: df.Function, eps: float, sign):\n",
" return (\n",
" params[\"mu1\"] * 0.5 * (1.0 + sign(l, eps))\n",
" + params[\"mu2\"] * 0.5 * (1.0 - sign(l, eps))\n",
" )"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"It is time to write down the equations.\n",
"\\begin{equation*}\n",
" \\mathbb{T} = \\mu (\\nabla v + (\\nabla v)^T) - p \\mathbb{I} \n",
"\\end{equation*}\n",
"\n",
"\\begin{equation*}\n",
" \\int_{\\Omega} \\rho \\partial_t v \\cdot \\varphi_{v} + \\rho ((\\nabla v) v) \\cdot \\varphi_v \\; dx + \\int_{\\Omega} \\mathbb{T} \\cdot \\nabla \\varphi_{v}\\; dx = 0\n",
"\\end{equation*}\n",
"\n",
"\\begin{equation*}\n",
" \\int_{\\Omega} \\text{div}(v)\\varphi_p \\; dx = 0\n",
"\\end{equation*}\n",
"\n",
"\\begin{equation*}\n",
" \\int_{\\Omega} (\\partial_t l) \\varphi_l + \\text{div}(l v) \\varphi_l \\;dx = 0 \n",
"\\end{equation*}"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"n = df.FacetNormal(mesh)\n",
"h = df.CellDiameter(mesh)\n",
"h_avg = (h('+') + h('-')) / 2.0\n",
"alpha = df.Constant(0.1)\n",
"\n",
"cauchy_green = (\n",
" mu(material_params, l, eps, sign) * (df.grad(v) + df.grad(v).T)\n",
" - p * df.Identity(mesh.topology().dim())\n",
")\n",
"\n",
"material_detivative = (\n",
" (1 / dt) * df.inner(v - v0, phi_v) # partial time derivative\n",
" + df.inner(df.grad(v) * v, phi_v) # convective therm\n",
")\n",
"\n",
"momentum = (\n",
" rho(material_params, l, eps, sign) * material_detivative*df.dx\n",
" + df.inner(cauchy_green, df.grad(phi_v))*df.dx\n",
" - rho(material_params, l, eps, sign)\n",
" * df.inner(g, phi_v) * df.dx\n",
")\n",
"\n",
"mass = (\n",
" df.div(v) * phi_p*df.dx\n",
")\n",
"\n",
"levelset_convection = (\n",
" (1 / dt) * df.inner(l - l0, phi_l) * df.dx\n",
" + df.div(l * v) * phi_l * df.dx\n",
")\n",
"\n",
"stabilization = (\n",
" alpha('+') * (h_avg**2)\n",
" * df.inner(df.jump(df.grad(l), n), df.jump(df.grad(phi_l), n)) * df.dS\n",
")\n",
"\n",
"pde_form = momentum + mass + levelset_convection + stabilization"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Our system is not linear, thus we need to solve it with Newton solver. Let's define it now!"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Set Newton-solver\n",
"# form compiler parameter\n",
"ffc_options = {\n",
" \"quadrature_degree\": 4,\n",
" \"optimize\": True,\n",
" \"eliminate_zeros\": True\n",
"}\n",
"\n",
"J = df.derivative(pde_form, w)\n",
"problem = df.NonlinearVariationalProblem(pde_form, w, bcs, J, ffc_options)\n",
"solver = df.NonlinearVariationalSolver(problem)\n",
"\n",
"prm = solver.parameters\n",
"prm['nonlinear_solver'] = 'newton'\n",
"prm['newton_solver']['linear_solver'] = 'mumps'\n",
"prm['newton_solver']['lu_solver']['report'] = False\n",
"prm['newton_solver']['absolute_tolerance'] = 1E-10\n",
"prm['newton_solver']['relative_tolerance'] = 1E-10\n",
"prm['newton_solver']['maximum_iterations'] = 20\n",
"prm['newton_solver']['report'] = True"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We create the XDMF files for storing the results."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# Initialize the files for writing the results\n",
"files = []\n",
"for name in ['v', 'p', 'l']:\n",
" with df.XDMFFile(df.MPI.comm_world, f\"result/{name}.xdmf\") as xdmf:\n",
" xdmf.parameters[\"flush_output\"] = True\n",
" files.append(xdmf)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"It remains to solve it!!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t = t_start\n",
"while t < t_end:\n",
" df.info(f\"t = {t}\")\n",
" solver.solve()\n",
" w0.assign(w)\n",
" t += dt\n",
" # write the time-step into the file\n",
" for func, name, xdmf in zip(w.split(True), ['v', 'p', 'l'], files):\n",
" func.rename(name, name)\n",
" xdmf.write(func, t)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Complete Code"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import dolfin as df\n",
"import matplotlib.pyplot as plt\n",
"\n",
"mesh = df.RectangleMesh(\n",
" df.Point(0.0, -0.5), df.Point(0.25, 0.5), 20, 80, 'crossed'\n",
")\n",
"df.plot(mesh)\n",
"plt.show()\n",
"\n",
"class InitialCondition(df.UserExpression):\n",
"\n",
" def __init__(self, **kwargs):\n",
" self.center = [0.125, 0.25]\n",
" self.radius = df.sqrt(0.125**2 + 0.25**2)\n",
" super().__init__(**kwargs)\n",
" \n",
" def r(self, x: list) -> float:\n",
" return df.sqrt((x[0] - self.center[0])**2 + (x[1] - self.center[1])**2)\n",
"\n",
" def eval(self, values, x):\n",
" values[0] = 0.0 # v_x\n",
" values[1] = 0.0 # v_y\n",
" values[2] = 0.0 # p\n",
" values[3] = self.r(x) - self.radius # l\n",
"\n",
" def value_shape(self):\n",
" return (4, )\n",
"# create instance of the class.\n",
"# initial_conditions = InitialCondition()\n",
"\n",
"\n",
"class InitialConditionSigmoid(df.UserExpression):\n",
"\n",
" def __init__(self, **kwargs):\n",
" self.center = [0.125, 0.25]\n",
" self.radius = df.sqrt(0.125**2 + 0.25**2)\n",
" self.eps = 0.001 \n",
" super().__init__(**kwargs)\n",
" \n",
" def r(self, x: list) -> float:\n",
" return df.sqrt((x[0] - self.center[0])**2 + (x[1] - self.center[1])**2)\n",
"\n",
" def eval(self, values, x):\n",
" values[0] = 0.0 # v_x\n",
" values[1] = 0.0 # v_y\n",
" values[2] = 0.0 # p\n",
" values[3] = 1 / (1 + df.exp(min((self.radius - self.r(x) ) / self.eps, 10)))\n",
" # values[3] = 1 / (df.exp( (self.r(x) - self.radius) / self.eps ))\n",
"\n",
" def value_shape(self):\n",
" return (4, )\n",
"\n",
"initial_conditions = InitialConditionSigmoid()\n",
"\n",
"# def sign(q: df.Function, eps: float):\n",
"# return q / df.sqrt(q * q + eps * eps * df.inner(df.grad(q), df.grad(q)))\n",
"\n",
"def sign(q: df.Function, eps: float):\n",
" return 2 * (q - 0.5) \n",
"\n",
"elements = [\n",
" df.VectorElement(\"CG\", mesh.ufl_cell(), 2), # velocity\n",
" df.FiniteElement(\"CG\", mesh.ufl_cell(), 1), # pressure\n",
" df.FiniteElement(\"CG\", mesh.ufl_cell(), 1) # levelset\n",
"]\n",
"\n",
"function_space = df.FunctionSpace(\n",
" mesh, df.MixedElement(elements)\n",
")\n",
"\n",
"material_params = {\n",
" \"mu1\": 1.0,\n",
" \"mu2\": 1.0,\n",
" \"rho1\": 500,\n",
" \"rho2\": 1000,\n",
"}\n",
"eps = 1e-4\n",
"dt = 0.02\n",
"t_start = 0.0\n",
"t_end = 1.0\n",
"g = df.Constant((0.0, -10.0)) # gravity field\n",
"\n",
"bcs = [\n",
" df.DirichletBC(\n",
" function_space.sub(0), # v\n",
" df.Constant((0.0, 0.0)),\n",
" \"on_boundary\"\n",
" )\n",
"]\n",
"\n",
"w = df.Function(function_space) # unknown\n",
"w0 = df.Function(function_space) # from previous step\n",
"phi = df.TestFunction(function_space) # test function\n",
"w0.assign(df.interpolate(initial_conditions, function_space))\n",
"w.assign(w0)\n",
"\n",
"# Split functions\n",
"v, p, l = df.split(w)\n",
"v0, p0, l0 = df.split(w0)\n",
"phi_v, phi_p, phi_l = df.split(phi)\n",
"\n",
"def rho(params: dict, l: df.Function, eps: float, sign):\n",
" return (\n",
" params[\"rho1\"] * 0.5* (1.0 + sign(l, eps))\n",
" + params[\"rho2\"] * 0.5 * (1.0 - sign(l, eps))\n",
" )\n",
"\n",
"\n",
"def mu(params: dict, l: df.Function, eps: float, sign):\n",
" return (\n",
" params[\"mu1\"] * 0.5 * (1.0 + sign(l, eps))\n",
" + params[\"mu2\"] * 0.5 * (1.0 - sign(l, eps))\n",
" )\n",
"\n",
"n = df.FacetNormal(mesh)\n",
"h = df.CellDiameter(mesh)\n",
"h_avg = (h('+') + h('-')) / 2.0\n",
"alpha = df.Constant(0.1)\n",
"\n",
"cauchy_green = (\n",
" mu(material_params, l, eps, sign) * (df.grad(v) + df.grad(v).T)\n",
" - p * df.Identity(mesh.topology().dim())\n",
")\n",
"\n",
"material_detivative = (\n",
" (1 / dt) * df.inner(v - v0, phi_v) # partial time derivative\n",
" + df.inner(df.grad(v) * v, phi_v) # convective therm\n",
")\n",
"\n",
"momentum = (\n",
" rho(material_params, l, eps, sign) * material_detivative*df.dx\n",
" + df.inner(cauchy_green, df.grad(phi_v))*df.dx\n",
" - rho(material_params, l, eps, sign)\n",
" * df.inner(g, phi_v) * df.dx\n",
")\n",
"\n",
"mass = (\n",
" df.div(v) * phi_p*df.dx\n",
")\n",
"\n",
"levelset_convection = (\n",
" (1 / dt) * df.inner(l - l0, phi_l) * df.dx\n",
" + df.div(l * v) * phi_l * df.dx\n",
")\n",
"\n",
"stabilization = (\n",
" alpha('+') * (h_avg**2)\n",
" * df.inner(df.jump(df.grad(l), n), df.jump(df.grad(phi_l), n)) * df.dS\n",
")\n",
"\n",
"pde_form = momentum + mass + levelset_convection + stabilization\n",
"\n",
"# Set Newton-solver\n",
"# form compiler parameter\n",
"ffc_options = {\n",
" \"quadrature_degree\": 4,\n",
" \"optimize\": True,\n",
" \"eliminate_zeros\": True\n",
"}\n",
"\n",
"J = df.derivative(pde_form, w)\n",
"problem = df.NonlinearVariationalProblem(pde_form, w, bcs, J, ffc_options)\n",
"solver = df.NonlinearVariationalSolver(problem)\n",
"\n",
"prm = solver.parameters\n",
"prm['nonlinear_solver'] = 'newton'\n",
"prm['newton_solver']['linear_solver'] = 'mumps'\n",
"prm['newton_solver']['lu_solver']['report'] = False\n",
"prm['newton_solver']['absolute_tolerance'] = 1E-10\n",
"prm['newton_solver']['relative_tolerance'] = 1E-10\n",
"prm['newton_solver']['maximum_iterations'] = 20\n",
"prm['newton_solver']['report'] = True\n",
"\n",
"# Initialize the files for writing the results\n",
"files = []\n",
"for name in ['v', 'p', 'l']:\n",
" with df.XDMFFile(df.MPI.comm_world, f\"result/{name}.xdmf\") as xdmf:\n",
" xdmf.parameters[\"flush_output\"] = True\n",
" files.append(xdmf)\n",
"\n",
"t = t_start\n",
"while t < t_end:\n",
" df.info(f\"t = {t}\")\n",
" solver.solve()\n",
" w0.assign(w)\n",
" t += dt\n",
" # write the time-step into the file\n",
" for func, name, xdmf in zip(w.split(True), ['v', 'p', 'l'], files):\n",
" func.rename(name, name)\n",
" xdmf.write(func, t)"
]
}
],
"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.10"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
![](\"pics/rt.png\")
\n",
"