Skip to content

Commit

Permalink
Merge pull request #152 from pmgbergen/robin_bc
Browse files Browse the repository at this point in the history 
Robin bc
  • Loading branch information
rbe051 authored Sep 14, 2018
2 parents 36b0d40 + 173beff commit 908d7a0
Show file tree
Hide file tree
Showing 23 changed files with 1,945 additions and 301 deletions.
184 changes: 180 additions & 4 deletions examples/example2/ConvergenceTestMPFA.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,10 +179,17 @@
"text": [
"Cartesian errors - no perturbations\n",
"[0.2177630338617197, 0.053854028644456195, 0.013418192291060563] pressure error\n",
"[4.04357927 4.01350849] convergence rate u^n-1 / u^n\n",
"\n",
"[0.020449627547754923, 0.00558959841388443, 0.0014594922987497814] flux error\n",
"[3.65851463 3.82982385] convergence rate f^n-1 / f^n\n",
"\n",
"Triangular errors - no perturbations\n",
"[0.1154963072658415, 0.028265747189434955, 0.007020068836548777] pressure error\n",
"[0.04223113257437278, 0.011882184577715638, 0.0033976662367782447] flux error\n"
"[4.08608718 4.02642023] convergence rate u^n-1 / u^n\n",
"\n",
"[0.04223113257437278, 0.011882184577715638, 0.0033976662367782447] flux error\n",
"[3.55415558 3.49716062] convergence rate f^n-1 / f^n\n"
]
}
],
Expand Down Expand Up @@ -276,15 +283,27 @@
"\n",
"\n",
"u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
"u_cart_rate = np.array(u_cart_nopert)\n",
"u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n",
"f_cart_rate = np.array(f_cart_nopert)\n",
"f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n",
"print('Cartesian errors - no perturbations')\n",
"print(u_cart_nopert,'pressure error')\n",
"print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n",
"print(f_cart_nopert,'flux error')\n",
"print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n",
"\n",
"u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
"u_triang_rate = np.array(u_triang_nopert)\n",
"u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n",
"f_triang_rate = np.array(f_triang_nopert)\n",
"f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n",
"\n",
"print('Triangular errors - no perturbations')\n",
"print(u_triang_nopert,'pressure error')\n",
"print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n",
"print(f_triang_nopert,'flux error')\n",
"\n",
"print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n",
"# EK: These values were hard-coded 11.06.2017.\n",
"assert np.abs(u_cart_nopert[2] - 0.01341819229106214) < 1e-10\n",
"assert np.abs(f_cart_nopert[2] - 0.0014594922987503009) < 1e-10\n",
Expand All @@ -294,6 +313,163 @@
"assert np.abs(f_triang_nopert[2] - 0.0033976662367782447) < 1e-10\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Robin boundary condition"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Cartesian errors - no perturbations\n",
"[0.562515011909662, 0.14015471392164033, 0.034984342689359696] pressure error\n",
"[4.01352902 4.00621258] convergence rate u^n-1 / u^n\n",
"\n",
"[0.22376060444435023, 0.061834410597352465, 0.01641158802403413] flux error\n",
"[3.61870684 3.76772866] convergence rate f^n-1 / f^n\n",
"\n",
"Triangular errors - no perturbations\n",
"[0.11192090940414492, 0.02516849817177693, 0.005995729167675974] pressure error\n",
"[4.44686483 4.19773767] convergence rate u^n-1 / u^n\n",
"\n",
"[0.03311849158106908, 0.00866805360137708, 0.002734187490448117] flux error\n",
"[3.82075297 3.17024843] convergence rate f^n-1 / f^n\n"
]
}
],
"source": [
"# Analytical solution\n",
"x, y = sympy.symbols('x y')\n",
"u = x * (1-x) * sympy.sin(x) * sympy.cos(y)\n",
"u_f = sympy.lambdify((x, y), u, 'numpy')\n",
"dux = sympy.diff(u, x)\n",
"duy = sympy.diff(u, y)\n",
"dux_f = sympy.lambdify((x, y), dux, 'numpy')\n",
"duy_f = sympy.lambdify((x, y), duy, 'numpy')\n",
"rhs = -sympy.diff(dux, x) - sympy.diff(duy, y)\n",
"rhs_f = sympy.lambdify((x, y), rhs, 'numpy')\n",
"\n",
"###\n",
"# This is where parameters can be modified to alter the convergence test.\n",
"# The remaining lines should \n",
"np.random.seed(42)\n",
"base = 4\n",
"domain = np.array([1, 1])\n",
"basedim = np.array([base, base])\n",
"num_refs = 3\n",
"grid_type = 'cart'\n",
"\n",
"def run_convergence(grid_type, pert):\n",
" u_err = []\n",
" flux_err = []\n",
"\n",
" for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
" # Reset the random seed for every grid realization.\n",
" # This should make no difference for the convergence test, \n",
" # but it makes sure that we can run unit tests based on the values obtained\n",
" # here.\n",
" np.random.seed(42)\n",
" \n",
" # Permeability tensor\n",
" k = perm(2, np.ones(g.num_cells))\n",
" robin_weight = 1.2\n",
" # Exact solution\n",
" xf = g.face_centers\n",
" xc = g.cell_centers\n",
" u_ex = u_f(xc[0], xc[1])\n",
" du_ex_faces = np.vstack((dux_f(xf[0], xf[1]), duy_f(xf[0], xf[1])))\n",
" flux_ex = -np.sum(g.face_normals[:2] * du_ex_faces, axis=0)\n",
"\n",
" # Set type of boundary conditions - Dirichlet\n",
" bound_faces = g.tags['domain_boundary_faces'].nonzero()[0]\n",
" n = g.nodes\n",
" top = bc.face_on_side(g, 'ymax')\n",
" bot = bc.face_on_side(g, 'ymin')\n",
" left = bc.face_on_side(g, 'xmin')\n",
" right = bc.face_on_side(g, 'xmax')\n",
"\n",
" dir_faces = np.asarray(left[0])\n",
" rob_faces = np.hstack((right[0]))\n",
" neu_faces = np.hstack((bot[0], top[0]))\n",
" bc_faces = np.hstack([dir_faces, rob_faces])\n",
" bc_names = ['dir']*len(dir_faces) + ['rob'] * len(rob_faces)\n",
" bound_cond = bc.BoundaryCondition(g, bc_faces, bc_names)\n",
" \n",
" # MPFA discretization, and system matrix\n",
" flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond, robin_weight=robin_weight)\n",
" div = fvutils.scalar_divergence(g)\n",
" a = div * flux\n",
" \n",
" # Boundary conditions\n",
" nfi, _, sgn_n = sps.find(g.cell_faces[neu_faces,:])\n",
" rfi, _, sgn = sps.find(g.cell_faces[rob_faces,:])\n",
" u_bound = np.zeros(g.num_faces)\n",
" u_bound[dir_faces] = u_f(xf[0, dir_faces], xf[1, dir_faces])\n",
" # innflow is always negative, so need to flip flux according to normal direction.\n",
" u_bound[neu_faces[nfi]] = flux_ex[neu_faces[nfi]] * (sgn_n)\n",
" u_bound[rob_faces[rfi]] = flux_ex[rob_faces[rfi]] * (sgn) +\\\n",
" robin_weight * u_f(xf[0, rob_faces], xf[1, rob_faces]) * g.face_areas[rob_faces]\n",
"\n",
" # Right hand side - contribution from the solution and the boundary conditions\n",
" xc = g.cell_centers\n",
" b = rhs_f(xc[0], xc[1]) * g.cell_volumes - div * bound_flux * u_bound\n",
"\n",
" # Solve system, derive fluxes\n",
" u_num = spsolve(a, b)\n",
" flux_num = flux * u_num + bound_flux * u_bound\n",
"\n",
" # Exact solution\n",
" flux_diff = flux_num - flux_ex\n",
" \n",
" u_err.append(np.sqrt(np.sum(g.cell_volumes * (u_num - u_ex)**2)) /\n",
" np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
" flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * flux_diff**2))/\n",
" np.sqrt(np.sum((g.face_areas ** g.dim) * flux_ex**2)))\n",
" return u_err, flux_err\n",
"\n",
"grids = ['cart', 'triangular']\n",
"\n",
"## No tests on perturbed grids - I'm too lazy to find faces of the respective boundaries\n",
"\n",
"\n",
"u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
"u_cart_rate = np.array(u_cart_nopert)\n",
"u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n",
"f_cart_rate = np.array(f_cart_nopert)\n",
"f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n",
"print('Cartesian errors - no perturbations')\n",
"print(u_cart_nopert,'pressure error')\n",
"print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n",
"print(f_cart_nopert,'flux error')\n",
"print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n",
"\n",
"u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
"u_triang_rate = np.array(u_triang_nopert)\n",
"u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n",
"f_triang_rate = np.array(f_triang_nopert)\n",
"f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n",
"\n",
"print('Triangular errors - no perturbations')\n",
"print(u_triang_nopert,'pressure error')\n",
"print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n",
"print(f_triang_nopert,'flux error')\n",
"print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n",
"\n",
"# RB: These values were hard-coded 11.09.2018.\n",
"assert np.abs(u_cart_nopert[2] - 0.034984342689359696) < 1e-10\n",
"assert np.abs(f_cart_nopert[2] - 0.01641158802403413) < 1e-10\n",
"assert np.abs(u_triang_nopert[2] - 0.005995729167675974) < 1e-10\n",
"assert np.abs(f_triang_nopert[2] - 0.002734187490448117) < 1e-10\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand All @@ -304,7 +480,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 5,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -471,7 +647,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 6,
"metadata": {},
"outputs": [
{
Expand Down
Loading

0 comments on commit 908d7a0

Please sign in to comment.