From 206cf8f44680266cdd2ef906ae95b2051c328fe1 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 7 Jun 2023 08:08:47 +0000 Subject: [PATCH 1/4] Fix warnings and parameter file for serial stokes. --- source/data_out.cc | 2 +- source/pdes/serial/stokes.cc | 11 +++++++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/source/data_out.cc b/source/data_out.cc index 6d5a2fe4a..890f1e2f1 100644 --- a/source/data_out.cc +++ b/source/data_out.cc @@ -177,7 +177,7 @@ namespace ParsedTools { static Vector material_ids; material_ids.reinit(dh.get_triangulation().n_active_cells()); - for (const auto cell : dh.active_cell_iterators()) + for (const auto &cell : dh.active_cell_iterators()) material_ids(cell->active_cell_index()) = cell->material_id(); data_out->add_data_vector( material_ids, diff --git a/source/pdes/serial/stokes.cc b/source/pdes/serial/stokes.cc index 657e7d38f..13e36df1c 100644 --- a/source/pdes/serial/stokes.cc +++ b/source/pdes/serial/stokes.cc @@ -52,7 +52,14 @@ namespace PDEs ParsedTools::Components::join(std::vector(dim + 1, "0"), ";"), "Exact solution") - , boundary_conditions("/Stokes/Boundary conditions", component_names) + , boundary_conditions( + "/Stokes/Boundary conditions", + component_names, + {{numbers::internal_face_boundary_id}}, + {"u"}, + {{ParsedTools::BoundaryConditionType::dirichlet}}, + {ParsedTools::Components::join(std::vector(dim + 1, "0"), + ";")}) , error_table(Utilities::split_string_list(component_names), {{VectorTools::H1_norm, VectorTools::L2_norm}, {VectorTools::L2_norm}}) @@ -60,7 +67,7 @@ namespace PDEs , velocity(0) , pressure(dim) { - enter_subsection("Error table"); + enter_subsection("Error table (" + component_names + ")"); enter_my_subsection(this->prm); error_table.add_parameters(this->prm); leave_my_subsection(this->prm); From a2c461035d6cc9175afc75c8446db5123d78bd20 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 7 Jun 2023 08:09:27 +0000 Subject: [PATCH 2/4] Added manufactured solution for Stokes. --- notebooks/manufactured_solution.ipynb | 132 ++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 notebooks/manufactured_solution.ipynb diff --git a/notebooks/manufactured_solution.ipynb b/notebooks/manufactured_solution.ipynb new file mode 100644 index 000000000..30086c7c7 --- /dev/null +++ b/notebooks/manufactured_solution.ipynb @@ -0,0 +1,132 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Manufactured solutions for Stokes problem\n", + "\n", + "When solving numerically the stokes problem \n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "-\\Delta \\boldsymbol{u} + \\nabla p & = f \\\\\n", + "\\nabla \\cdot \\boldsymbol{u} &= 0,\n", + "\\end{aligned}\n", + "$$\n", + "\n", + "we can construct manufactured solutions that satisfy the divergence free\n", + "constraint, by taking a solution $\\boldsymbol{u}$ that is the curl of a scalar\n", + "field (in two-dimensions) or the curl of a vector field in three dimensions.\n", + "\n", + "For two dimensions: \n", + "\n", + "$$\n", + "\\begin{split}\n", + "\\boldsymbol{u}_x & = \\partial_y g \\\\\n", + "\\boldsymbol{u}_y & = -\\partial_x g\n", + "\\end{split}\n", + "$$\n", + "\n", + "which satisfy by construction $\\nabla \\cdot \\boldsymbol{u} = 0$ for $g\\in\n", + "C^2(\\Omega)$.\n", + "\n", + "We can then build a forcing term $f$ that would force the system to have this\n", + "exact solution, by using the equality $f = -\\Delta \\boldsymbol{u} + \\nabla p$.\n", + "\n", + "If we want to use these in input files for the `stokes` executable of the\n", + "`fsi-suite`, we must replace the `**` with `^`, and we must write them in a\n", + "format that deal.II understands. This is done in the following cells, where we\n", + "first compute symbolically the differential operators, and then produce an\n", + "output that is readable by deal.II." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import *\n", + "import textwrap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def write_parameters(precursor, p):\n", + " # ux, uy, and p are the exact solution of the velocity and pressure\n", + " ux = diff(precursor, y)\n", + " uy = -diff(precursor, x)\n", + "\n", + " # lap_u_x and lap_u_y are minus the laplacian of ux and uy\n", + " lap_u_x = - diff(ux, x, 2) - diff(ux, y, 2)\n", + " lap_u_y = - diff(uy, x, 2) - diff(uy, y, 2)\n", + "\n", + " # fx and fy are the resulting forcing terms\n", + " fx = lap_u_x + diff(p, x)\n", + " fy = lap_u_y + diff(p, y)\n", + "\n", + " def prm(ux, uy, p, fx, fy):\n", + " def to_prm(x):\n", + " text = str(x).replace('**', '^')\n", + " start_len = len(\" set Forcing term = \")\n", + " return ' \\\\\\n'.join(textwrap.wrap(text, 80-start_len, break_long_words=False, \n", + " subsequent_indent=' '*start_len))\n", + "\n", + " print(\n", + " \"subsection Functions\\n set Exact solution = \", to_prm(ux), \";\",\n", + " to_prm(uy), \";\", to_prm(p), \"\\n\",\n", + " \" set Forcing term = \", to_prm(fx), \";\",\n", + " to_prm(fy), \"; 0\", \"\\nend\\n\"\n", + " )\n", + " \n", + " prm(ux, uy, p, fx, fy)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Spatial variables\n", + "x,y = symbols('x y')\n", + "\n", + "# Precursor of th exact solution. We take the curl of this function to have a \n", + "# divergence free velocity field\n", + "precursor = sin(pi*x)**2*sin(pi*y)**2\n", + "pressure = cos(2*pi*x)*cos(2*pi*y)\n", + "\n", + "write_parameters(precursor, pressure)" + ] + } + ], + "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.15" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 707748ab389afb5967a424da323a01a6ae30b6c8 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 10 Jul 2023 22:20:05 +0200 Subject: [PATCH 3/4] Single brace on scalars. --- source/pdes/linear_visco_elasticity.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/pdes/linear_visco_elasticity.cc b/source/pdes/linear_visco_elasticity.cc index 9c04224da..223289686 100644 --- a/source/pdes/linear_visco_elasticity.cc +++ b/source/pdes/linear_visco_elasticity.cc @@ -44,7 +44,7 @@ namespace PDEs "Second Lame coefficient", "Shear viscosity", "Bulk viscosity"}) - , material_ids_0({{0}}) + , material_ids_0({0}) , eulerian_mapping(this->dof_handler, "/LinearViscoElasticity/Mapping") { this->add_parameter("Material ids of region 0", material_ids_0); From a14f66f0a6b83157ad784ce53c3d453334cf303c Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 25 Jun 2024 02:21:34 +0200 Subject: [PATCH 4/4] logstream. --- include/parsed_tools/grid_generator.h | 1 + source/grid_generator.cc | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/include/parsed_tools/grid_generator.h b/include/parsed_tools/grid_generator.h index 262340360..ceef85596 100644 --- a/include/parsed_tools/grid_generator.h +++ b/include/parsed_tools/grid_generator.h @@ -16,6 +16,7 @@ #ifndef parsed_tools_grid_generator_h #define parsed_tools_grid_generator_h +#include #include #include diff --git a/source/grid_generator.cc b/source/grid_generator.cc index a0d547915..089d0e544 100644 --- a/source/grid_generator.cc +++ b/source/grid_generator.cc @@ -275,4 +275,4 @@ namespace ParsedTools template class GridGenerator<2, 3>; template class GridGenerator<3, 3>; -} // namespace ParsedTools \ No newline at end of file +} // namespace ParsedTools