Earthquake_Deformation_v1.0.1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "55dd5a17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# SRC comes from this link: https://fenicsproject.org/pub/tutorial/html/._ftut1008.html\n",
    "# Version 1.0.0\n",
    "# 11/23/2021: version 1.0.0. The code can simulate deformation subjuected to far-fielding shear loading. (D.Liu)\n",
    "\n",
    "# Elastic deformation from 3D faulting simulated by modulus reduction in the fault zone.\n",
    "# Created on 11/23/2021 by Dunyu Liu (dliu@ig.utexas.edu)\n",
    "# Modified some parameters 11/27/2021 by Xu (xiaohua.xu@austin.utexas.edu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fe01794a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.\n",
      "WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.\n",
      "Assembling the system and preconditioner ...\n",
      "Solving linear system of size 509343 x 509343 (PETSc Krylov solver).\n",
      "Solving the system and writing out the results ...\n",
      "  0 KSP preconditioned resid norm 2.533300796643e+01 true resid norm 1.047765230940e+14 ||r(i)||/||b|| 1.000000000000e+00\n",
      "  1 KSP preconditioned resid norm 3.209449355397e+00 true resid norm 8.137921758567e+12 ||r(i)||/||b|| 7.766932437017e-02\n",
      "  2 KSP preconditioned resid norm 7.533841740519e-01 true resid norm 3.276072884913e+12 ||r(i)||/||b|| 3.126724182262e-02\n",
      "  3 KSP preconditioned resid norm 9.440307262528e-02 true resid norm 4.523459509904e+11 ||r(i)||/||b|| 4.317245291530e-03\n",
      "  4 KSP preconditioned resid norm 1.372092079581e-02 true resid norm 5.290254665612e+10 ||r(i)||/||b|| 5.049084002211e-04\n",
      "  5 KSP preconditioned resid norm 4.306900875256e-03 true resid norm 9.817677772594e+09 ||r(i)||/||b|| 9.370112199454e-05\n",
      "  6 KSP preconditioned resid norm 1.270006638917e-03 true resid norm 4.279207363016e+09 ||r(i)||/||b|| 4.084128043815e-05\n",
      "  7 KSP preconditioned resid norm 2.412656658742e-04 true resid norm 6.937641634085e+08 ||r(i)||/||b|| 6.621370350170e-06\n",
      "  8 KSP preconditioned resid norm 3.724126162240e-05 true resid norm 1.162931185394e+08 ||r(i)||/||b|| 1.109915800843e-06\n",
      "  9 KSP preconditioned resid norm 6.790448084056e-06 true resid norm 2.286472311146e+07 ||r(i)||/||b|| 2.182237245165e-07\n",
      " 10 KSP preconditioned resid norm 1.173262143517e-06 true resid norm 2.851446681556e+06 ||r(i)||/||b|| 2.721455720570e-08\n",
      " 11 KSP preconditioned resid norm 3.029279784426e-07 true resid norm 7.890296649315e+05 ||r(i)||/||b|| 7.530595992612e-09\n",
      " 12 KSP preconditioned resid norm 4.534345064039e-08 true resid norm 1.335773726813e+05 ||r(i)||/||b|| 1.274878844390e-09\n",
      " 13 KSP preconditioned resid norm 7.130400662266e-09 true resid norm 2.243342842017e+04 ||r(i)||/||b|| 2.141073950321e-10\n",
      "\n",
      "UMFPACK V5.7.9 (Oct 20, 2019): ERROR: out of memory\n",
      "\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "\n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to successfully call PETSc function 'KSPSolve'.\n*** Reason:  PETSc error code is: 76 (Error in external library).\n*** Where:   This error was encountered inside /Users/runner/miniforge3/conda-bld/fenics-pkgs_1637474655529/work/dolfin/dolfin/la/PETScKrylovSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2019.1.0\n*** Git changeset:  \n*** -------------------------------------------------------------------------\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/0w/shxyptc10k1b2zyhgh5prrnw0000gq/T/ipykernel_9528/1264575905.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    198\u001b[0m \u001b[0mvon_Mises\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3.\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0minner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    199\u001b[0m \u001b[0mV\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFunctionSpace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'P'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 200\u001b[0;31m \u001b[0mvon_Mises\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mproject\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvon_Mises\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mV\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    201\u001b[0m \u001b[0;31m#plot(von_Mises, title='Stress intensity')\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    202\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/homebrew/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/fem/projection.py\u001b[0m in \u001b[0;36mproject\u001b[0;34m(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)\u001b[0m\n\u001b[1;32m    136\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mfunction\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    137\u001b[0m         \u001b[0mfunction\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mV\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 138\u001b[0;31m     \u001b[0mcpp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mla\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvector\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msolver_type\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpreconditioner_type\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    139\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    140\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: \n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to successfully call PETSc function 'KSPSolve'.\n*** Reason:  PETSc error code is: 76 (Error in external library).\n*** Where:   This error was encountered inside /Users/runner/miniforge3/conda-bld/fenics-pkgs_1637474655529/work/dolfin/dolfin/la/PETScKrylovSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2019.1.0\n*** Git changeset:  \n*** -------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Linear elastic problem.\n",
    "\n",
    "  -div(sigma(u)) = f\n",
    "\n",
    "The model is used to simulate earthquake deformation. \n",
    "Faulting is simulated via modulus reduction in the fault zone. \n",
    "\"\"\"\n",
    "\n",
    "from __future__ import print_function\n",
    "from fenics import *\n",
    "from ufl import nabla_div\n",
    "\n",
    "# Model domain.\n",
    "L = 5e3; W = 5e3; Depth = 5e3\n",
    "n1, n2, n3 = 100, 40, 40\n",
    "\n",
    "# Rock properties\n",
    "#vp = 6000\n",
    "#vs = 3464\n",
    "rho = 2670\n",
    "#mu = vs*vs*rho\n",
    "mu = 20e9\n",
    "nu = 0.25\n",
    "# shear_modulus inside the fault zone\n",
    "mu_fault = mu*0.5\n",
    "\n",
    "lambda_ = 2*mu*nu/(1-2*nu)\n",
    "lambda_fault = 2*mu_fault*nu/(1-2*nu)\n",
    "g = 0 # gravitiy, assumed to be zero in this application.\n",
    "\n",
    "\n",
    "shear_slip = 0.05 # far-field loading by slip\n",
    "normal_slip = -0.1\n",
    "fault_zone_width = 1e3 # fault zone width\n",
    "\n",
    "dim = 3\n",
    "\n",
    "# Create mesh and define function space\n",
    "mesh = BoxMesh(Point(0, 0, 0), Point(L, W, Depth), n1, n2, n3)\n",
    "V = VectorFunctionSpace(mesh, 'P', 1)\n",
    "boundaries = MeshFunction(\"size_t\", mesh, dim-1)\n",
    "\n",
    "tol = 1E-14\n",
    "\n",
    "# Define boundary condition\n",
    "class Top_boundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return near(x[2], Depth, DOLFIN_EPS) and on_boundary\n",
    "class Bottom_boundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return near(x[2], 0.0, DOLFIN_EPS) and on_boundary\n",
    "class Left_boundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return near(x[0], 0, DOLFIN_EPS) and on_boundary\n",
    "class Right_boundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):   \n",
    "        return near(x[0], L, DOLFIN_EPS) and on_boundary\n",
    "class North_boundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return near(x[1], W, DOLFIN_EPS) and on_boundary\n",
    "class South_boundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return near(x[1], 0, DOLFIN_EPS) and on_boundary    \n",
    "boundaries.set_all(0)\n",
    "Top_boundary().mark(boundaries, 1)\n",
    "Bottom_boundary().mark(boundaries, 2)\n",
    "Left_boundary().mark(boundaries, 3)\n",
    "Right_boundary().mark(boundaries, 4)\n",
    "North_boundary().mark(boundaries, 5)\n",
    "South_boundary().mark(boundaries, 6)\n",
    "# Rename boundaries\n",
    "top = 1\n",
    "bottom = 2\n",
    "left = 3\n",
    "right = 4\n",
    "north = 5\n",
    "south = 6\n",
    "# setting boundary conditions\n",
    "#bc = [DirichletBC(V.sub(2), Constant((0.0)), boundaries, bottom),\n",
    "#        DirichletBC(V.sub(0), Constant((0.0)), boundaries, left),\n",
    "#        DirichletBC(V.sub(0), Constant((0.0)), boundaries, right),\n",
    "#        DirichletBC(V.sub(1), shear_slip, boundaries, left), # far-field loading via slip\n",
    "#        DirichletBC(V.sub(1), -shear_slip, boundaries, right)] # far-field loading via slip\n",
    "\n",
    "#bc = [  DirichletBC(V.sub(0), Constant((0.0)), boundaries, left),\n",
    "#        DirichletBC(V.sub(0), Constant((0.0)), boundaries, right),\n",
    "#        DirichletBC(V.sub(1), shear_slip, boundaries, left), # far-field loading via slip\n",
    "#        DirichletBC(V.sub(1), -shear_slip, boundaries, right)] # far-field loading via slip\n",
    "\n",
    "bc = [  DirichletBC(V.sub(0), -normal_slip, boundaries, left),\n",
    "        DirichletBC(V.sub(0), normal_slip, boundaries, right),\n",
    "        DirichletBC(V.sub(2), Constant((0.0)), boundaries, bottom)] \n",
    "\n",
    "#bc = [DirichletBC(V.sub(2), Constant((0.0)), boundaries, bottom),\n",
    "#        DirichletBC(V.sub(0), Constant((0.0)), boundaries, left),\n",
    "#        DirichletBC(V.sub(0), normal_slip, boundaries, right),\n",
    "#        DirichletBC(V.sub(1), Constant((0.0)), boundaries, left), # far-field loading via slip\n",
    "#        DirichletBC(V.sub(1), Constant((0.0)), boundaries, right)] # far-field loading via slip\n",
    "\n",
    "\n",
    "# Define fault zone as subdomains.\n",
    "# At the center of x axis.\n",
    "class Omega_0(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return x[0] >= L/2-fault_zone_width/2 and x[0] <= L/2+fault_zone_width/2 \\\n",
    "            and x[1]> W/2-W/3 and x[1] < W/2+W/3 and  x[2] >= fault_zone_width\n",
    "mf = MeshFunction(\"size_t\", mesh, 3)\n",
    "subdomain0 = Omega_0()\n",
    "#subdomain1 = Omega_1()\n",
    "mf.set_all(10)\n",
    "subdomain0.mark(mf,11) # 1, the middle anisotrpic layer.\n",
    "rock = 10\n",
    "fault = 11\n",
    "\n",
    "# Define strain and stress\n",
    "def epsilon(u):\n",
    "    return 0.5*(nabla_grad(u) + nabla_grad(u).T)\n",
    "    #return sym(nabla_grad(u))\n",
    "\n",
    "def sigma(u):\n",
    "    return lambda_mod*nabla_div(u)*Identity(d) + 2*shear_mod*epsilon(u)\n",
    "\n",
    "# Define class to set different parameters to different subdomains\n",
    "class K(UserExpression):\n",
    "            def __init__(self, subdomains, k_fault, k_rock, **kwargs):\n",
    "                super().__init__(**kwargs)\n",
    "                self.subdomains = subdomains\n",
    "                self.k_fault = k_fault\n",
    "                self.k_rock = k_rock\n",
    "\n",
    "            def eval_cell(self, values, x, cell):\n",
    "                if self.subdomains[cell.index] == fault: #anisotrpic layer\n",
    "                    values[0] = self.k_fault\n",
    "                elif self.subdomains[cell.index] == rock:\n",
    "                    values[0] = self.k_rock\n",
    "            def values_shape(self):\n",
    "                return (1,)            \n",
    "shear_mod = K(mf, mu_fault, mu, degree=0)\n",
    "lambda_mod = K(mf,lambda_fault, lambda_, degree=0)\n",
    "\n",
    "# Define variational problem\n",
    "u = TrialFunction(V)\n",
    "d = u.geometric_dimension()  # space dimension\n",
    "v = TestFunction(V)\n",
    "f = Constant((0, 0, -rho*g))\n",
    "T = Constant((0, 0, 0))\n",
    "a = inner(sigma(u), epsilon(v))*dx\n",
    "b = inner(sigma(u), epsilon(v))*dx\n",
    "L = dot(f, v)*dx + dot(T, v)*ds\n",
    "\n",
    "# Compute solution\n",
    "#u = Function(V)\n",
    "#solve(a == L, u, bc)\n",
    "\n",
    "A, bb = assemble_system(a, L, bc)\n",
    "#A, bb = assemble_system(a, L)\n",
    "#for bc in bcs:\n",
    "#    bc.apply(A)\n",
    "#    bc.apply(bb)\n",
    "# Assemble preconditioner system\n",
    "P, btmp = assemble_system(b, L, bc)\n",
    "#P, btmp = assemble_system(b, L)\n",
    "#for bc in bcs:\n",
    "#    bc.apply(P)\n",
    "print('Assembling the system and preconditioner ...')\n",
    " \n",
    "# Create Krylov solver and AMG preconditioner\n",
    " \n",
    " \n",
    "#solver = KrylovSolver(\"minres\", \"amg\") # best.\n",
    "solver = KrylovSolver(\"gmres\", \"amg\") # second best.\n",
    "#solver = KrylovSolver(\"tfqmr\", \"amg\") # least\n",
    " \n",
    "solver.parameters[\"relative_tolerance\"] = 1.0e-9   # most important\n",
    "solver.parameters[\"absolute_tolerance\"] = 1.0e-15  # set a small number\n",
    "solver.parameters[\"divergence_limit\"] = 1.0e4\n",
    "solver.parameters[\"maximum_iterations\"] = 7000\n",
    "solver.parameters[\"error_on_nonconvergence\"] = True\n",
    "solver.parameters[\"nonzero_initial_guess\"] = False\n",
    "solver.parameters[\"report\"] = True\n",
    "solver.parameters[\"monitor_convergence\"] = True\n",
    "# Associate operator (A) and preconditioner matrix (P)\n",
    "#solver.set_operators(A, P)\n",
    "solver.set_operators(A, P)\n",
    "\n",
    "# Solve\n",
    "u = Function(V)\n",
    "solver.solve(u.vector(), bb)\n",
    "print('Solving the system and writing out the results ...')\n",
    "\n",
    "\n",
    "# Plot solution\n",
    "#plot(u, title='Displacement', mode='displacement')\n",
    "\n",
    "# Plot stress\n",
    "s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress\n",
    "von_Mises = sqrt(3./2*inner(s, s))\n",
    "V = FunctionSpace(mesh, 'P', 1)\n",
    "von_Mises = project(von_Mises, V)\n",
    "#plot(von_Mises, title='Stress intensity')\n",
    "\n",
    "# Compute magnitude of displacement\n",
    "u_magnitude = sqrt(dot(u, u))\n",
    "u_magnitude = project(u_magnitude, V)\n",
    "#plot(u_magnitude, 'Displacement magnitude')\n",
    "#print('min/max u:',\n",
    "#      u_magnitude.vector().array().min(),\n",
    "#      u_magnitude.vector().array().max())\n",
    "\n",
    "# Save solution to file in VTK format\n",
    "File('elasticity/displacement.pvd') << u\n",
    "#File('elasticity/von_mises.pvd') << von_Mises\n",
    "#File('elasticity/magnitude.pvd') << u_magnitude\n",
    "\n",
    "# Hold plot\n",
    "#interactive()\n",
    "print('Simulation Done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aee5566a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "fenicsproject",
   "language": "python",
   "name": "fenicsproject"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}