RS IN JULIA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "94a41dea",
   "metadata": {},
   "source": [
    "# Method of Regularized Stokeslet\n",
    "## Introduction\n",
    "**This work was done as a part of the final project for a special topics course titled _Boundary integral method_ in the fall of 22 by [Professor Cortez](https://sse.tulane.edu/math/faculty/cortez)**\n",
    "\n",
    "The fluid flow problems in tiny scales are usually modeled  by the Stokes equations for incompressible flows \n",
    "  $$\\mu \\Delta \\textbf{u}=\\nabla p-\\textbf{F}$$\n",
    "  $$ \\nabla \\cdot \\textbf{u}=0$$\n",
    "where $\\mu$ is the fluid viscosity, $p$ is the pressure, $\\textbf{u}$ is the velocity, and $\\textbf{F}$ is force. A fundamental solution of these equations is called a $Stokeslet$. The particular case of a single force $\\mathbf{f}_0$ exerted at $\\mathbf{x}_0$ results in a velocity  field given by\n",
    "$$\\mathbf{u}=\\frac{\\mathbf{f}_0}{8\\pi\\mu r}+\\frac{(\\mathbf{f_0\\cdot(x-x_0)})(\\mathbf{x-x_0})}{8\n",
    "\\pi\\mu r}$$\n",
    "where $r=||\\mathbf{x-x_0}||$.\n",
    "\n",
    "Note this solution is undefined at $r=0$ or $\\mathbf{x=x}_0$.\n",
    "\n",
    "However, the singularities can be eliminated through the function(usually known as blob function) $\\phi_{\\delta}(\\mathbf{x})$ which is radially symmetric and satisfies that the integral over the space is one. So, considering $\\mathbf{F}=\\mathbf{f}_0\\phi_{\\delta}$ the singularity can be removed. The idea is due to [Professor Cortez](https://epubs.siam.org/doi/10.1137/S106482750038146X).\n",
    "\n",
    "With the following choice of blob function:\n",
    "$$\\phi_\\delta(r)=\\dfrac{15\\delta^4}{8\\pi(r^2+\\delta^2)^{\\frac{7}{2}}}$$\n",
    "the regularized stokeslet is \n",
    "$$\\vec{u}(x)=\\boldsymbol f_0 \\dfrac{r^2+2\\delta^2}{8\\pi\\mu(r^2+\\delta^2)^{\\frac{3}{2}}}+\\dfrac{(\\boldsymbol f_{0} .\\boldsymbol x)\\boldsymbol x}{8\\pi\\mu(r^2+\\delta^2)^{\\frac{3}{2}}}$$\n",
    "\n",
    "\n",
    "\n",
    "## Swimming of filament immersed in a viscous fluid\n",
    "\n",
    "We suppose the slender body is a sine wave,\n",
    "$$y(s)=A\\cos(\\lambda s-2\\pi t),z(s)=0$$\n",
    "and $x(s)$ such that $\\sqrt{(x')^2+(y')^2}=1$,\n",
    "and the curvature \n",
    "\\begin{align*}\n",
    "    \\kappa=&\\frac{x''y'-y''x'}{(\\sqrt{(x')^2+(y')^2})^{3}}\\\\=&x''y'-y''x'\\\\=&\\frac{-y''}{\\sqrt{1-(y')^2}}\n",
    "\\end{align*}\n",
    "To compute the forces, we've utilized the approach discussed in the article [A computational model of aquatic animal locomotion\n",
    "](https://www.sciencedirect.com/science/article/abs/pii/0021999188901581).\n",
    "\n",
    "Once we've computed forces at each point, we can calculate velocity. Then, the locomotion of filament can be found by solving $\\dfrac{d\\boldsymbol X}{dt}=\\boldsymbol u$. To solve this we have used the forward Euler method(other methods like RK can also be used).\n",
    "\n",
    "It's worth mentioning that this model has been utilized for various problems; especially those related to the [motility of hyper-activated mammalian sperm](https://www.sciencedirect.com/science/article/abs/pii/S0022519314001635). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "28be53ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "velocityRS (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "function velocityRS(x, y, s, F, delta, δs)\n",
    "\n",
    "    u = zeros(2, length(s))\n",
    "    int_u = zeros(size(u))\n",
    "    r1 = zeros(size(x))\n",
    "    r2 = zeros(size(y))\n",
    "    F_dot_r = zeros(size(x))\n",
    "    R2 = zeros(size(x))\n",
    "    for i = 1:length(s)\n",
    "\n",
    "        r1 .= x .- x[i]\n",
    "        r2 .= y .- y[i]\n",
    "        R2 .= r1 .^ 2 + r2 .^ 2\n",
    "        F_dot_r .= (F[1, i] .* r1) .+ (F[2, i] .* r2)\n",
    "\n",
    "        int_u[1, :] = (1 / (8 * pi)) * ((F[1, i] .* (R2 .+ 2 * delta^2)) ./ ((R2 .+ delta^2) .^ 1.5) .+ (F_dot_r .* (r1)) ./ ((R2 .+ delta^2) .^ 1.5))\n",
    "\n",
    "        int_u[2, :] = (1 / (8 * pi)) * ((F[2, i] .* (R2 .+ 2 * delta^2)) ./ ((R2 .+ delta^2) .^ 1.5) .+ (F_dot_r .* (r2)) ./ ((R2 .+ delta^2) .^ 1.5))\n",
    "\n",
    "        u[1, :] = u[1, :] .+ (int_u[1, :])\n",
    "        u[2, :] = u[2, :] .+ (int_u[2, :])\n",
    "\n",
    "    end\n",
    "    return u\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1bb16ceb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Forces (generic function with 2 methods)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function Forces(x,y,s,δs,t0,k,S_T,S_C)\n",
    "    F = zeros(2, length(s));\n",
    "\n",
    "    for i = 1:length(s)-1\n",
    "        X_s = [x[i+1] .- x[i], y[i+1] .- y[i]]\n",
    "        X_s_distance2 = (X_s[1] ^ 2 + X_s[2] ^ 2)\n",
    "        F_t = -S_T * ((sqrt(X_s_distance2) / δs) - 1) * (X_s / (sqrt(X_s_distance2)))\n",
    "\n",
    "\n",
    "        F[1, i] = F[1, i] - F_t[1]\n",
    "        F[1, i+1] = F[1, i+1] + F_t[1]\n",
    "        F[2, i] = F[2, i] - F_t[2]\n",
    "        F[2, i+1] = F[2, i+1] + F_t[2]\n",
    "    end\n",
    "\n",
    "    for i = 2:length(s)-1\n",
    "        κ = k.(s[i], t0)\n",
    "\n",
    "        D2_x_i = (x[i+1] .- 2 * x[i] .+ x[i-1]) / (δs^2)\n",
    "        D0_x_i = (x[i+1] .- x[i-1]) / (2 * δs)\n",
    "\n",
    "        D2_y_i = (y[i+1] .- 2 * y[i] .+ y[i-1]) / (δs^2)\n",
    "        D0_y_i = (y[i+1] .- y[i-1]) / (2 * δs)\n",
    "\n",
    "\n",
    "        F1 = S_C * ((D2_x_i  * D0_y_i  - D0_x_i  * D2_y_i) - κ) * (D0_y_i / δs)\n",
    "        F2 = -S_C * ((D2_x_i * D0_y_i - D0_x_i * D2_y_i) - κ) * (D0_x_i / δs)\n",
    "        Q1 = -S_C * 0.5 * ((D2_x_i * D0_y_i - D0_x_i * D2_y_i) - κ) * (D2_y_i)\n",
    "        Q2 = S_C * 0.5 * ((D2_x_i * D0_y_i - D0_x_i * D2_y_i) - κ) * (D2_x_i)\n",
    "\n",
    "\n",
    "\n",
    "        F[1, i+1] = F[1, i+1] - F1 - Q1\n",
    "        F[1, i] = F[1, i] + 2 * F1\n",
    "        F[1, i-1] = F[1, i-1] - F1 + Q1\n",
    "\n",
    "        F[2, i+1] = F[2, i+1] - F2 - Q2\n",
    "        F[2, i] = F[2, i] .+ 2 * F2\n",
    "        F[2, i-1] = F[2, i-1] - F2 + Q2\n",
    "    end\n",
    "    return F\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "69c5691a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RSFILL (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# import Pkg; Pkg.add(\"DifferentialEquations\")\n",
    "# import Pkg; Pkg.add(\"Plots\")\n",
    "using DifferentialEquations, Plots\n",
    "# params = (S_T = 75,S_C = 0.075,L = 7,Amp = 0.1,N = 51, dt = 1e-4)\n",
    "\n",
    "function RSFILL(S_T, S_C, L, Amp, N, t0, tmax, dt)\n",
    "\n",
    "\n",
    "    λ = 2 * π / L\n",
    "    δs = L / (N - 1)\n",
    "    s_range = (0.0, L)\n",
    "\n",
    "    δ =  δs\n",
    "    f(u, p, s) = sqrt(1 - Amp^2 * λ^2 * (cos(λ * s)) .^ 2)\n",
    "\n",
    "    x0 = 0\n",
    "    prob = ODEProblem(f, x0, s_range)\n",
    "    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat=δs)\n",
    "\n",
    "    y = Amp * sin.(λ * sol.t)\n",
    "\n",
    "\n",
    "    delta = (2 / 3) * δs\n",
    "\n",
    "    s = LinRange(0, L, N)\n",
    "    x = sol'\n",
    "    y = y'\n",
    "\n",
    "    k(ss, tt) = Amp * λ^2 * sin(λ * ss - 2 * pi * tt) / (sqrt(1 - Amp^2 * λ^2 * (cos(λ * ss - 2 * pi * tt))^2))\n",
    "\n",
    "    nmax = round((tmax - t0) / dt)\n",
    "    anim = Animation()\n",
    "    for n = 1:nmax\n",
    "\n",
    "        F = Forces(x, y, s, δs, t0, k, S_T, S_C)\n",
    "\n",
    "        u = velocityRS(x, y, s, F, delta, δs)\n",
    "\n",
    "        for i = 1:length(s)\n",
    "\n",
    "            x[i] = x[i] + dt * u[1, i]\n",
    "            y[i] = y[i] + dt * u[2, i]\n",
    "\n",
    "        end\n",
    "\n",
    "        t0 = t0 + dt\n",
    "\n",
    "        if mod(n, 100) == 0\n",
    "\n",
    "            plot(x', y', xlims=(-L, L), ylims=(-2*Amp, 2*Amp), aspect_ratio=:equal, legend=:false, marker=:circle, lw=3)\n",
    "            #     quiver!(x',y',quiver=(F[1,:]',F[2,:]'))\n",
    "\n",
    "            annotate!(0.5, 0.25, \"Time= $(rpad(round(t0; digits=2),4,\"0\")) s\")\n",
    "            frame(anim)\n",
    "\n",
    "        end\n",
    "\n",
    "    end\n",
    "\n",
    "    vid = gif(anim, \"RS1.gif\", fps=10)\n",
    "    return vid\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f454636c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Saved animation to c:\\Users\\mosle\\RegularizedStokeslet\\RS1.gif\n",
      "└ @ Plots C:\\Users\\mosle\\.julia\\packages\\Plots\\sxUvK\\src\\animation.jl:156\n"
     ]
    },
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: anim not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: anim not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\Users\\mosle\\RegularizedStokeslet\\RS IN JULIA.ipynb:3"
     ]
    }
   ],
   "source": [
    "# RSFILL(S_T, S_C, L, Amp, N, t0, tmax, dt)\n",
    "RSFILL(75,  0.05, 1,  0.1, 31,0,3, 1e-5)\n"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": "dfbbe9037b684dff8d4e06f3c48e206c",
   "lastKernelId": "7ec23082-e3be-43f5-8e30-2f975259ee7a"
  },
  "kernelspec": {
   "display_name": "Julia 1.8.5",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}