bi_saturation_main.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from SpinBiSaturation.BiSaturation import BiSaturation\n",
    "import numpy as np\n",
    "from bvpdae_solver import BVPDAE, BVPSol\n",
    "from timeit import default_timer as dtime\n",
    "from utils_functions import plot_sequence\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eps =  0.1\n",
      "alpha =  0.5\n",
      "     # Residual error = 1.5873525112487923 with N = 501\n",
      "     # Residual error = 7.962064672998207 with N = 635\n",
      "     # Residual error = 250.48166379579965 with N = 1517\n",
      "     # Residual error = 587.3089346827488 with N = 3816\n",
      "     # Residual error = 2517.086341151253 with N = 9035\n",
      "     # Residual error = 2530.1202551662727 with N = 20881\n",
      "     # Residual error = 3158.3242768943396 with N = 49029\n",
      "Solving complete\n",
      "infos NLSE =  1.7202308334469312\n",
      "iter =  0\n"
     ]
    },
    {
     "ename": "Exception",
     "evalue": "Primal-dual BiSaturation problem failed before convergence",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mException\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[14], line 84\u001b[0m\n\u001b[0;32m     80\u001b[0m     plt\u001b[38;5;241m.\u001b[39mshow()\n\u001b[0;32m     83\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;18m__name__\u001b[39m \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m__main__\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m---> 84\u001b[0m     \u001b[43mmain_primal_dual\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[1;32mIn[14], line 26\u001b[0m, in \u001b[0;36mmain_primal_dual\u001b[1;34m()\u001b[0m\n\u001b[0;32m     24\u001b[0m cumul_time \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m ts1 \u001b[38;5;241m-\u001b[39m ts0\n\u001b[0;32m     25\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m infos\u001b[38;5;241m.\u001b[39msuccess:\n\u001b[1;32m---> 26\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPrimal-dual BiSaturation problem failed before convergence\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m     27\u001b[0m time, xp, z \u001b[38;5;241m=\u001b[39m bvp_sol\u001b[38;5;241m.\u001b[39mtime, bvp_sol\u001b[38;5;241m.\u001b[39mxp, bvp_sol\u001b[38;5;241m.\u001b[39mz\n\u001b[0;32m     28\u001b[0m convergence \u001b[38;5;241m=\u001b[39m eps \u001b[38;5;241m<\u001b[39m\u001b[38;5;241m=\u001b[39m tol\n",
      "\u001b[1;31mException\u001b[0m: Primal-dual BiSaturation problem failed before convergence"
     ]
    }
   ],
   "source": [
    "def main_primal_dual():\n",
    "    cumul_time = 0.\n",
    "    eps = .1\n",
    "    alpha = 0.5\n",
    "    tol = 1e-7\n",
    "    ocp = BiSaturation()\n",
    "    ocp.set_eps(eps)\n",
    "    options = dict(control_odes_error=True, display=1, res_tol=1e-3)\n",
    "    bvp_solver = BVPDAE(**options)\n",
    "    time, xp, z = ocp.initialize()\n",
    "    initial_solution = dict(time=time, xp=xp, z=z)\n",
    "    convergence = False\n",
    "    iter = 0\n",
    "    times, x0s, x1s, x2s, x3s, p0s, p1s, p2s, p3s, us = [list() for _ in range(10)]\n",
    "    bvp_sol = BVPSol(time=time, xp=xp, z=z)\n",
    "    while not convergence:\n",
    "        print(\"eps = \", eps)\n",
    "        print(\"alpha = \", alpha)\n",
    "        ts0 = dtime()\n",
    "        bvp_sol, infos = bvp_solver.solve(bvp_sol, ocp)\n",
    "        print(\"infos NLSE = \", infos.NLSE_infos.bc_residual)\n",
    "        print(\"iter = \", iter)\n",
    "        ts1 = dtime()\n",
    "        cumul_time += ts1 - ts0\n",
    "        if not infos.success:\n",
    "            raise Exception(\"Primal-dual BiSaturation problem failed before convergence\")\n",
    "        time, xp, z = bvp_sol.time, bvp_sol.xp, bvp_sol.z\n",
    "        convergence = eps <= tol\n",
    "        eps *= alpha\n",
    "        ocp.set_eps(eps)\n",
    "        times.append(xp[4, 0] * time)\n",
    "        x0s.append(xp[0, :])\n",
    "        x1s.append(xp[1, :])\n",
    "        x2s.append(xp[2, :])\n",
    "        x3s.append(xp[3, :])\n",
    "        p0s.append(xp[5, :])\n",
    "        p1s.append(xp[6, :])\n",
    "        p2s.append(xp[7, :])\n",
    "        p3s.append(xp[8, :])\n",
    "        us.append(z[0])\n",
    "        iter += 1\n",
    "        print(\" \")\n",
    "    optimal_solution = dict(time=time, xp=xp, z=z, eps=eps / alpha, alpha=alpha, exec_time=cumul_time, iter=iter)\n",
    "    dict_save = dict(initial_solution=initial_solution, optimal_solution=optimal_solution)\n",
    "    with open(\"SpinBiSaturation/results_BiSaturation_primal_dual.pickle\", \"wb\") as fh:\n",
    "        pickle.dump(dict_save, fh)\n",
    "    plot_sequence(times, x0s, \"$x_0(t)$\",\n",
    "                  \"Sequence of optimal penalized $\\\\bar{x}_0$\")\n",
    "    plot_sequence(times, x1s, \"$x_1(t)$\",\n",
    "                  \"Sequence of optimal penalized $\\\\bar{x}_1$\")\n",
    "    plot_sequence(times, x2s, \"$x_2(t)$\",\n",
    "                    \"Sequence of optimal penalized $\\\\bar{x}_2$\")\n",
    "    plot_sequence(times, x3s, \"$x_3(t)$\",\n",
    "                    \"Sequence of optimal penalized $\\\\bar{x}_3$\")\n",
    "    plot_sequence(times, p0s, \"$p_0(t)$\",\n",
    "                  \"Sequence of optimal penalized $\\\\bar{p}_0$\")\n",
    "    plot_sequence(times, p1s, \"$p_1(t)$\",\n",
    "                  \"Sequence of optimal penalized $\\\\bar{p}_1$\")\n",
    "    plot_sequence(times, p2s, \"$p_2(t)$\",\n",
    "                    \"Sequence of optimal penalized $\\\\bar{p}_2$\")\n",
    "    plot_sequence(times, p3s, \"$p_3(t)$\",\n",
    "                    \"Sequence of optimal penalized $\\\\bar{p}_3$\")\n",
    "    plot_sequence(times, us, \"$u(t)$\",\n",
    "                  \"Sequence of optimal penalized $\\\\bar{u}_0$\")\n",
    " \n",
    "\n",
    "    t = np.linspace(0., 2. * np.pi, 1000)\n",
    "    x1t = ocp.r * ocp.a1 * np.cos(t) + ocp.xfinal[0] / 2.\n",
    "    x2t = ocp.r * ocp.a2 * np.sin(t) + ocp.xfinal[1] / 2.5\n",
    "    plt.figure()\n",
    "    plt.plot(xp[0], xp[1])\n",
    "    plt.plot(x1t, x2t)\n",
    "\n",
    "    plt.figure()\n",
    "\n",
    "    print(\"Time final = \", xp[4, 0])\n",
    "    print(\"Execution time = \", cumul_time)\n",
    "    print(\"Number of iterations = \", len(times))\n",
    "    print(\"Number of final time steps = \", len(time))\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main_primal_dual()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "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.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}