From 3175f8251058fa6e7686810e2320947bddc6ff8c Mon Sep 17 00:00:00 2001 From: David Zwicker Date: Wed, 5 Feb 2025 10:04:40 +0100 Subject: [PATCH] Fixed 9-point stencil in Cartesian coordinates --- docs/methods/stencils/Cartesian2D.nb | 564 +++++++++++++++++++++++++++ pde/grids/boundaries/axis.py | 2 +- pde/grids/operators/cartesian.py | 37 +- pde/tools/config.py | 4 +- tests/grids/test_cartesian_grids.py | 14 +- 5 files changed, 597 insertions(+), 24 deletions(-) create mode 100644 docs/methods/stencils/Cartesian2D.nb diff --git a/docs/methods/stencils/Cartesian2D.nb b/docs/methods/stencils/Cartesian2D.nb new file mode 100644 index 00000000..d59f1708 --- /dev/null +++ b/docs/methods/stencils/Cartesian2D.nb @@ -0,0 +1,564 @@ +(* Content-type: application/vnd.wolfram.mathematica *) + +(*** Wolfram Notebook File ***) +(* http://www.wolfram.com/nb *) + +(* CreatedBy='Wolfram 14.1' *) + +(*CacheID: 234*) +(* Internal cache information: +NotebookFileLineBreakTest +NotebookFileLineBreakTest +NotebookDataPosition[ 154, 7] +NotebookDataLength[ 19313, 556] +NotebookOptionsPosition[ 16928, 509] +NotebookOutlinePosition[ 17327, 525] +CellTagsIndexPosition[ 17284, 522] +WindowFrame->Normal*) + +(* Beginning of Notebook Content *) +Notebook[{ +Cell[BoxData[ + RowBox[{"Import", "[", + "\"\\"", "]"}]], "Input", + CellChangeTimes->{{3.947331747822672*^9, 3.947331747824991*^9}}, + CellLabel->"In[47]:=",ExpressionUUID->"02ae3299-3cef-44e9-bd9d-26279c613410"], + +Cell[BoxData[ + RowBox[{ + RowBox[{"stencilIso", "[", + RowBox[{"\[Gamma]_", ",", + RowBox[{"\[Delta]_", ":", "1"}]}], "]"}], ":=", + RowBox[{ + RowBox[{"(", + RowBox[{ + RowBox[{ + RowBox[{"(", + RowBox[{"1", "-", "\[Gamma]"}], ")"}], "*", + RowBox[{"{", + RowBox[{ + RowBox[{"{", + RowBox[{"0", ",", "1", ",", "0"}], "}"}], ",", + RowBox[{"{", + RowBox[{"1", ",", + RowBox[{"-", "4"}], ",", "1"}], "}"}], ",", + RowBox[{"{", + RowBox[{"0", ",", "1", ",", "0"}], "}"}]}], "}"}]}], "+", + RowBox[{"\[Gamma]", "*", + RowBox[{"{", + RowBox[{ + RowBox[{"{", + RowBox[{ + RowBox[{"1", "/", "2"}], ",", "0", ",", + RowBox[{"1", "/", "2"}]}], "}"}], ",", + RowBox[{"{", + RowBox[{"0", ",", + RowBox[{"-", "2"}], ",", "0"}], "}"}], ",", + RowBox[{"{", + RowBox[{ + RowBox[{"1", "/", "2"}], ",", "0", ",", + RowBox[{"1", "/", "2"}]}], "}"}]}], "}"}]}]}], ")"}], "/", + RowBox[{"\[Delta]", "^", "2"}]}]}]], "Input", + CellChangeTimes->{{3.947331332313459*^9, 3.947331393328499*^9}, { + 3.94733156280276*^9, 3.947331585253502*^9}}, + CellLabel->"In[48]:=",ExpressionUUID->"bbca8448-63fb-481c-acd5-85d0bad438bf"], + +Cell[BoxData[ + RowBox[{ + RowBox[{"stencilGen", "[", + RowBox[{"\[Gamma]_", ",", + RowBox[{"\[Delta]x_", ":", "1"}], ",", + RowBox[{"\[Delta]y_", ":", "1"}]}], "]"}], ":=", + RowBox[{"Module", "[", + RowBox[{ + RowBox[{"{", "stencil1D", "}"}], ",", "\[IndentingNewLine]", + RowBox[{ + RowBox[{"stencil1D", "=", + RowBox[{ + RowBox[{ + RowBox[{"(", + RowBox[{"1", "-", "\[Gamma]"}], ")"}], "*", + RowBox[{"{", + RowBox[{ + RowBox[{"{", + RowBox[{"0", ",", "1", ",", "0"}], "}"}], ",", + RowBox[{"{", + RowBox[{"0", ",", + RowBox[{"-", "2"}], ",", "0"}], "}"}], ",", + RowBox[{"{", + RowBox[{"0", ",", "1", ",", "0"}], "}"}]}], "}"}]}], "+", + RowBox[{"\[Gamma]", "*", + RowBox[{"{", + RowBox[{ + RowBox[{"{", + RowBox[{ + RowBox[{"1", "/", "4"}], ",", "0", ",", + RowBox[{"1", "/", "4"}]}], "}"}], ",", + RowBox[{"{", + RowBox[{"0", ",", + RowBox[{"-", "1"}], ",", "0"}], "}"}], ",", + RowBox[{"{", + RowBox[{ + RowBox[{"1", "/", "4"}], ",", "0", ",", + RowBox[{"1", "/", "4"}]}], "}"}]}], "}"}]}]}]}], ";", + "\[IndentingNewLine]", + RowBox[{ + RowBox[{"stencil1D", "/", + RowBox[{"\[Delta]x", "^", "2"}]}], "+", + RowBox[{ + RowBox[{"Transpose", "[", "stencil1D", "]"}], "/", + RowBox[{"\[Delta]y", "^", "2"}]}]}]}]}], "\[IndentingNewLine]", + "]"}]}]], "Input", + CellChangeTimes->{{3.947331395175108*^9, 3.947331526740199*^9}, { + 3.9473315882422028`*^9, 3.947331622544409*^9}, {3.9475804680440273`*^9, + 3.947580476816614*^9}, {3.947580725372374*^9, 3.947580734611315*^9}}, + CellLabel->"In[49]:=",ExpressionUUID->"746552fa-8f11-46be-9498-9060eb1e8ef6"], + +Cell[BoxData[ + RowBox[{ + RowBox[{"stencilGen2", "[", + RowBox[{"\[Gamma]_", ",", + RowBox[{"\[Delta]x_", ":", "1"}], ",", + RowBox[{"\[Delta]y_", ":", "1"}]}], "]"}], ":=", + RowBox[{"Module", "[", + RowBox[{ + RowBox[{"{", + RowBox[{"stencil", "=", + RowBox[{"ConstantArray", "[", + RowBox[{"0", ",", + RowBox[{"{", + RowBox[{"3", ",", "3"}], "}"}]}], "]"}]}], "}"}], ",", + "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"2", ",", "2"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"(", + RowBox[{"\[Gamma]", "-", "2"}], ")"}], "*", + RowBox[{"(", + RowBox[{ + RowBox[{"\[Delta]x", "^", + RowBox[{"-", "2"}]}], "+", + RowBox[{"\[Delta]y", "^", + RowBox[{"-", "2"}]}]}], ")"}]}]}], ";", "\[IndentingNewLine]", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"1", ",", "2"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"3", ",", "2"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"(", + RowBox[{"1", "-", "\[Gamma]"}], ")"}], "/", + RowBox[{"\[Delta]x", "^", "2"}]}]}]}], ";", "\[IndentingNewLine]", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"2", ",", "1"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"2", ",", "3"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"(", + RowBox[{"1", "-", "\[Gamma]"}], ")"}], "/", + RowBox[{"\[Delta]y", "^", "2"}]}]}]}], ";", "\[IndentingNewLine]", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"1", ",", "1"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"1", ",", "3"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"3", ",", "3"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{"stencil", "[", + RowBox[{"[", + RowBox[{"3", ",", "1"}], "]"}], "]"}], "=", + RowBox[{ + RowBox[{ + RowBox[{"\[Gamma]", "/", "\[Delta]x"}], "/", "\[Delta]y"}], "/", + "2"}]}]}]}]}], ";", "\[IndentingNewLine]", "stencil"}]}], + "\[IndentingNewLine]", "]"}]}]], "Input", + CellChangeTimes->{{3.947586908678421*^9, 3.947586925139866*^9}, { + 3.94758695855357*^9, 3.9475870565087433`*^9}, {3.947587102708673*^9, + 3.9475871562743*^9}, {3.947587340761941*^9, 3.947587358287702*^9}}, + CellLabel->"In[88]:=",ExpressionUUID->"68a058b4-4af3-46cc-ad2d-5d5afe23e012"], + +Cell[BoxData[ + RowBox[{ + RowBox[{"stencil5", "[", + RowBox[{"\[Delta]x_", ",", "\[Delta]y_"}], "]"}], ":=", + RowBox[{"{", + RowBox[{ + RowBox[{"{", + RowBox[{"0", ",", + RowBox[{"\[Delta]x", "^", + RowBox[{"-", "2"}]}], ",", "0"}], "}"}], ",", + RowBox[{"{", + RowBox[{ + RowBox[{"\[Delta]y", "^", + RowBox[{"-", "2"}]}], ",", + RowBox[{ + RowBox[{"-", "2"}], "*", + RowBox[{"(", + RowBox[{ + RowBox[{"\[Delta]x", "^", + RowBox[{"-", "2"}]}], "+", + RowBox[{"\[Delta]y", "^", + RowBox[{"-", "2"}]}]}], ")"}]}], ",", + RowBox[{"\[Delta]y", "^", + RowBox[{"-", "2"}]}]}], "}"}], ",", + RowBox[{"{", + RowBox[{"0", ",", + RowBox[{"\[Delta]x", "^", + RowBox[{"-", "2"}]}], ",", "0"}], "}"}]}], "}"}]}]], "Input", + CellChangeTimes->{{3.947568775425036*^9, 3.947568840882539*^9}, { + 3.9475804801741056`*^9, 3.947580487074244*^9}}, + CellLabel->"In[89]:=",ExpressionUUID->"fc0eec24-cbf2-4650-a212-c3fa881a96c2"], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"stencilIso", "[", + RowBox[{"\[Gamma]", ",", "\[Delta]"}], "]"}], "==", + RowBox[{"stencilGen2", "[", + RowBox[{"\[Gamma]", ",", "\[Delta]", ",", "\[Delta]"}], "]"}]}], + "]"}]], "Input", + CellChangeTimes->{{3.947331485848175*^9, 3.947331494037428*^9}, { + 3.947331531242619*^9, 3.9473315400423307`*^9}, {3.947331617104992*^9, + 3.947331651567712*^9}, 3.9475871613296413`*^9}, + CellLabel->"In[90]:=",ExpressionUUID->"08715573-1df8-42ca-ab53-06c273ea1af8"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{ + 3.947331494321125*^9, {3.9473315334561863`*^9, 3.947331540497634*^9}, { + 3.947331628651463*^9, 3.947331656869886*^9}, 3.9473317514380913`*^9, + 3.9475804889900637`*^9, 3.947580736306797*^9, 3.947587161993099*^9, + 3.9475873984026527`*^9}, + CellLabel->"Out[90]=",ExpressionUUID->"5be01652-f8db-4d0c-bd8d-6dda431f5c1c"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"FullSimplify", "[", + RowBox[{ + RowBox[{"stencilGen2", "[", + RowBox[{"0", ",", "\[Delta]x", ",", "\[Delta]y"}], "]"}], "==", + RowBox[{"stencil5", "[", + RowBox[{"\[Delta]x", ",", "\[Delta]y"}], "]"}]}], "]"}]], "Input", + CellChangeTimes->{{3.947568844037901*^9, 3.947568876807642*^9}, + 3.94758716506501*^9}, + CellLabel->"In[91]:=",ExpressionUUID->"69126aa7-899f-4912-b405-e9e1c5591403"], + +Cell[BoxData["True"], "Output", + CellChangeTimes->{{3.9475688659275723`*^9, 3.947568877063353*^9}, + 3.947580489386592*^9, 3.947580736316513*^9, 3.947587165659683*^9, + 3.947587398778366*^9}, + CellLabel->"Out[91]=",ExpressionUUID->"4a696e72-9766-42a5-9274-841ca32abbfd"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{ + RowBox[{ + RowBox[{"stencilGen2", "[", + RowBox[{"0", ",", "dx", ",", "dy"}], "]"}], "//", "FullSimplify"}], "//", + "MatrixForm"}]], "Input", + CellChangeTimes->{3.947580916143017*^9, 3.9475871675848093`*^9}, + CellLabel->"In[92]:=",ExpressionUUID->"7fc88a9f-5d01-44c8-b04f-90b462b4ed1c"], + +Cell[BoxData[ + TagBox[ + RowBox[{"(", "\[NoBreak]", GridBox[{ + {"0", + FractionBox["1", + SuperscriptBox["dx", "2"]], "0"}, + { + FractionBox["1", + SuperscriptBox["dy", "2"]], + RowBox[{ + RowBox[{"-", "2"}], " ", + RowBox[{"(", + RowBox[{ + FractionBox["1", + SuperscriptBox["dx", "2"]], "+", + FractionBox["1", + SuperscriptBox["dy", "2"]]}], ")"}]}], + FractionBox["1", + SuperscriptBox["dy", "2"]]}, + {"0", + FractionBox["1", + SuperscriptBox["dx", "2"]], "0"} + }, + GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}}, + GridBoxSpacings->{"Columns" -> { + Offset[0.27999999999999997`], { + Offset[0.7]}, + Offset[0.27999999999999997`]}, "Rows" -> { + Offset[0.2], { + Offset[0.4]}, + Offset[0.2]}}], "\[NoBreak]", ")"}], + Function[BoxForm`e$, + MatrixForm[BoxForm`e$]]]], "Output", + CellChangeTimes->{3.947580916378909*^9, 3.947587167908883*^9, + 3.947587399206608*^9}, + CellLabel-> + "Out[92]//MatrixForm=",ExpressionUUID->"492ea3c6-d74f-4876-bdab-\ +cf3e8eb7cd38"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{ + RowBox[{ + RowBox[{"stencilGen2", "[", + RowBox[{"w", ",", "dx", ",", "dx"}], "]"}], "//", "FullSimplify"}], "//", + "MatrixForm"}]], "Input", + CellChangeTimes->{{3.947586950290687*^9, 3.947586950392149*^9}, + 3.947587176020821*^9}, + CellLabel->"In[93]:=",ExpressionUUID->"03596475-eb3e-42a0-8832-37616b5d81e7"], + +Cell[BoxData[ + TagBox[ + RowBox[{"(", "\[NoBreak]", GridBox[{ + { + FractionBox["w", + RowBox[{"2", " ", + SuperscriptBox["dx", "2"]}]], + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dx", "2"]], + FractionBox["w", + RowBox[{"2", " ", + SuperscriptBox["dx", "2"]}]]}, + { + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dx", "2"]], + FractionBox[ + RowBox[{"2", " ", + RowBox[{"(", + RowBox[{ + RowBox[{"-", "2"}], "+", "w"}], ")"}]}], + SuperscriptBox["dx", "2"]], + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dx", "2"]]}, + { + FractionBox["w", + RowBox[{"2", " ", + SuperscriptBox["dx", "2"]}]], + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dx", "2"]], + FractionBox["w", + RowBox[{"2", " ", + SuperscriptBox["dx", "2"]}]]} + }, + GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}}, + GridBoxSpacings->{"Columns" -> { + Offset[0.27999999999999997`], { + Offset[0.7]}, + Offset[0.27999999999999997`]}, "Rows" -> { + Offset[0.2], { + Offset[0.4]}, + Offset[0.2]}}], "\[NoBreak]", ")"}], + Function[BoxForm`e$, + MatrixForm[BoxForm`e$]]]], "Output", + CellChangeTimes->{3.947586950679983*^9, 3.9475871765990343`*^9, + 3.947587399875432*^9}, + CellLabel-> + "Out[93]//MatrixForm=",ExpressionUUID->"3b7ea4d8-ec40-4463-8e88-\ +c7f6e7561052"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{ + RowBox[{ + RowBox[{"stencilGen2", "[", + RowBox[{"w", ",", "dx", ",", "dy"}], "]"}], "//", "FullSimplify"}], "//", + "MatrixForm"}]], "Input", + CellChangeTimes->{{3.947580753508809*^9, 3.947580804750575*^9}, + 3.947587180378683*^9}, + CellLabel->"In[94]:=",ExpressionUUID->"bee98d28-8ac8-4190-bd36-c0d2bcab2d3c"], + +Cell[BoxData[ + TagBox[ + RowBox[{"(", "\[NoBreak]", GridBox[{ + { + FractionBox["w", + RowBox[{"2", " ", "dx", " ", "dy"}]], + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dx", "2"]], + FractionBox["w", + RowBox[{"2", " ", "dx", " ", "dy"}]]}, + { + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dy", "2"]], + RowBox[{ + RowBox[{"(", + RowBox[{ + FractionBox["1", + SuperscriptBox["dx", "2"]], "+", + FractionBox["1", + SuperscriptBox["dy", "2"]]}], ")"}], " ", + RowBox[{"(", + RowBox[{ + RowBox[{"-", "2"}], "+", "w"}], ")"}]}], + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dy", "2"]]}, + { + FractionBox["w", + RowBox[{"2", " ", "dx", " ", "dy"}]], + FractionBox[ + RowBox[{"1", "-", "w"}], + SuperscriptBox["dx", "2"]], + FractionBox["w", + RowBox[{"2", " ", "dx", " ", "dy"}]]} + }, + GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}}, + GridBoxSpacings->{"Columns" -> { + Offset[0.27999999999999997`], { + Offset[0.7]}, + Offset[0.27999999999999997`]}, "Rows" -> { + Offset[0.2], { + Offset[0.4]}, + Offset[0.2]}}], "\[NoBreak]", ")"}], + Function[BoxForm`e$, + MatrixForm[BoxForm`e$]]]], "Output", + CellChangeTimes->{{3.947580755818836*^9, 3.947580807462764*^9}, + 3.947587180919794*^9, 3.9475874002356586`*^9}, + CellLabel-> + "Out[94]//MatrixForm=",ExpressionUUID->"73e99b65-aa2d-4843-9e1d-\ +c0f5bf89de52"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"FullSimplify", "@", + RowBox[{"Total", "@", + RowBox[{"Flatten", "@", + RowBox[{"stencilGen2", "[", + RowBox[{"w", ",", "dx", ",", "dy"}], "]"}]}]}]}]], "Input", + CellChangeTimes->{{3.94758729636064*^9, 3.947587316866667*^9}}, + CellLabel->"In[79]:=",ExpressionUUID->"4e498f19-06f5-4de2-8a2b-e4ae7b9865f5"], + +Cell[BoxData[ + RowBox[{"-", + FractionBox[ + RowBox[{ + SuperscriptBox[ + RowBox[{"(", + RowBox[{"dx", "-", "dy"}], ")"}], "2"], " ", "w"}], + RowBox[{ + SuperscriptBox["dx", "2"], " ", + SuperscriptBox["dy", "2"]}]]}]], "Output", + CellChangeTimes->{{3.947587304772172*^9, 3.947587317112425*^9}}, + CellLabel->"Out[79]=",ExpressionUUID->"d8e2f0f4-9aa5-4396-95b9-aee912f6b294"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"ToPython", "[", + RowBox[{ + RowBox[{ + RowBox[{"stencilGen", "[", + RowBox[{"w", ",", "dx", ",", "dy"}], "]"}], "//", "FullSimplify"}], ",", + + RowBox[{"Copy", "->", "True"}]}], "]"}]], "Input", + CellChangeTimes->{{3.9473316594642897`*^9, 3.947331673233478*^9}, { + 3.947331754688772*^9, 3.947331811476253*^9}}, + CellLabel->"In[53]:=",ExpressionUUID->"cb21bde0-8d42-4687-8416-7755020f6946"], + +Cell[BoxData["\<\"np.array([np.array([0.25 * (dx ** -2 + dy ** -2) * w, (dx \ +** -2) * (1 -w), 0.25 * (dx ** -2 + dy ** -2) * w]), np.array([(dy ** -2) * \ +(1 -w), (dx ** -2) * (dy ** -2) * (dx ** 2 + dy ** 2) * (-2 + w), (dy ** -2) \ +* (1 -w)]), np.array([0.25 * (dx ** -2 + dy ** -2) * w, (dx ** -2) * (1 -w), \ +0.25 * (dx ** -2 + dy ** -2) * w])])\"\>"], "Output", + CellChangeTimes->{{3.9473316663492002`*^9, 3.94733167342844*^9}, + 3.9473317159945793`*^9, {3.9473317566700563`*^9, 3.9473318118993683`*^9}, + 3.947580490014529*^9, 3.947580736336179*^9}, + CellLabel->"Out[53]=",ExpressionUUID->"d8282a14-7050-49d4-9526-bceb671a25e0"] +}, Open ]], + +Cell[BoxData[""], "Input", + CellChangeTimes->{{3.947580491246051*^9, 3.9475804917531652`*^9}}, + CellLabel->"In[54]:=",ExpressionUUID->"d278b8c0-bcc1-4bf2-a18f-494250c612f6"] +}, +WindowSize->{1512, 916}, +WindowMargins->{{-1512, Automatic}, {-271, Automatic}}, +FrontEndVersion->"14.1 for Mac OS X ARM (64-bit) (July 16, 2024)", +StyleDefinitions->"Default.nb", +ExpressionUUID->"9194eac2-888c-4fff-a307-e37036fe28b1" +] +(* End of Notebook Content *) + +(* Internal cache information *) +(*CellTagsOutline +CellTagsIndex->{} +*) +(*CellTagsIndex +CellTagsIndex->{} +*) +(*NotebookFileOutline +Notebook[{ +Cell[554, 20, 302, 5, 29, "Input",ExpressionUUID->"02ae3299-3cef-44e9-bd9d-26279c613410"], +Cell[859, 27, 1289, 37, 29, "Input",ExpressionUUID->"bbca8448-63fb-481c-acd5-85d0bad438bf"], +Cell[2151, 66, 1847, 49, 90, "Input",ExpressionUUID->"746552fa-8f11-46be-9498-9060eb1e8ef6"], +Cell[4001, 117, 2725, 77, 151, "Input",ExpressionUUID->"68a058b4-4af3-46cc-ad2d-5d5afe23e012"], +Cell[6729, 196, 1030, 30, 29, "Input",ExpressionUUID->"fc0eec24-cbf2-4650-a212-c3fa881a96c2"], +Cell[CellGroupData[{ +Cell[7784, 230, 535, 11, 29, "Input",ExpressionUUID->"08715573-1df8-42ca-ab53-06c273ea1af8"], +Cell[8322, 243, 377, 6, 33, "Output",ExpressionUUID->"5be01652-f8db-4d0c-bd8d-6dda431f5c1c"] +}, Open ]], +Cell[CellGroupData[{ +Cell[8736, 254, 426, 9, 29, "Input",ExpressionUUID->"69126aa7-899f-4912-b405-e9e1c5591403"], +Cell[9165, 265, 274, 4, 33, "Output",ExpressionUUID->"4a696e72-9766-42a5-9274-841ca32abbfd"] +}, Open ]], +Cell[CellGroupData[{ +Cell[9476, 274, 318, 7, 29, "Input",ExpressionUUID->"7fc88a9f-5d01-44c8-b04f-90b462b4ed1c"], +Cell[9797, 283, 1148, 37, 104, "Output",ExpressionUUID->"492ea3c6-d74f-4876-bdab-cf3e8eb7cd38"] +}, Open ]], +Cell[CellGroupData[{ +Cell[10982, 325, 344, 8, 29, "Input",ExpressionUUID->"03596475-eb3e-42a0-8832-37616b5d81e7"], +Cell[11329, 335, 1537, 51, 101, "Output",ExpressionUUID->"3b7ea4d8-ec40-4463-8e88-c7f6e7561052"] +}, Open ]], +Cell[CellGroupData[{ +Cell[12903, 391, 344, 8, 29, "Input",ExpressionUUID->"bee98d28-8ac8-4190-bd36-c0d2bcab2d3c"], +Cell[13250, 401, 1606, 51, 107, "Output",ExpressionUUID->"73e99b65-aa2d-4843-9e1d-c0f5bf89de52"] +}, Open ]], +Cell[CellGroupData[{ +Cell[14893, 457, 339, 7, 29, "Input",ExpressionUUID->"4e498f19-06f5-4de2-8a2b-e4ae7b9865f5"], +Cell[15235, 466, 394, 11, 52, "Output",ExpressionUUID->"d8e2f0f4-9aa5-4396-95b9-aee912f6b294"] +}, Open ]], +Cell[CellGroupData[{ +Cell[15666, 482, 428, 10, 29, "Input",ExpressionUUID->"cb21bde0-8d42-4687-8416-7755020f6946"], +Cell[16097, 494, 639, 8, 53, "Output",ExpressionUUID->"d8282a14-7050-49d4-9526-bceb671a25e0"] +}, Open ]], +Cell[16751, 505, 173, 2, 29, "Input",ExpressionUUID->"d278b8c0-bcc1-4bf2-a18f-494250c612f6"] +} +] +*) + diff --git a/pde/grids/boundaries/axis.py b/pde/grids/boundaries/axis.py index 73b24bc8..d1640040 100644 --- a/pde/grids/boundaries/axis.py +++ b/pde/grids/boundaries/axis.py @@ -303,7 +303,7 @@ def from_data(cls, grid: GridBase, axis: int, data, rank: int = 0) -> BoundaryPa except TypeError: # if len is not supported, the format must be wrong raise BCDataError( - f"Unsupported boundary format: `{data}`." + cls.get_help() + f"Unsupported boundary format: `{data}`. " + cls.get_help() ) from None else: if data_len == 2: diff --git a/pde/grids/operators/cartesian.py b/pde/grids/operators/cartesian.py index de7dc7ec..7dc1d49a 100644 --- a/pde/grids/operators/cartesian.py +++ b/pde/grids/operators/cartesian.py @@ -362,8 +362,8 @@ def _make_laplace_numba_2d( grid (:class:`~pde.grids.cartesian.CartesianGrid`): The grid for which the operator is created corner_weight (float): - Weighting factor for the corner points. If `None`, the value is read from - the configuration option `operators.cartesian_laplacian_2d_corner_weight`. + Weighting factor for the corner points of the stencil. If `None`, the value + is read from the configuration option `operators.cartesian_laplacian_2d_corner_weight`. The standard value is zero, which corresponds to the traditional 5-point stencil. Typical alternative choices are 1/2 (Oono-Puri stencil) and 1/3 (Patra-Karttunen or Mehrstellen stencil); see @@ -380,7 +380,7 @@ def _make_laplace_numba_2d( parallel = dim_x * dim_y >= config["numba.multithreading_threshold"] if corner_weight == 0: - # 5-point stencil + # use standard 5-point stencil scale_x, scale_y = grid.discretization**-2 @jit(parallel=parallel) @@ -393,28 +393,25 @@ def laplace(arr: np.ndarray, out: np.ndarray) -> None: out[i - 1, j - 1] = lap_x + lap_y else: - # 9-point stencil + # use 9-point stencil with interpolated boundary conditions w = corner_weight _logger.info("Create 2D Cartesian Laplacian with 9-point stencil (w=%.3g)", w) - dx, dy = grid.discretization**-2 + if not np.isclose(*grid.discretization): + # we have not yet found a good expression for the 9-point stencil for dx!=dy + _logger.warning( + "9-point stencils with anisotropic grids are not tested and might " + "produce wrong results." + ) + + # prepare the stencil matrix + dxm2, dym2 = grid.discretization**-2 + dm2 = dxm2 + dym2 stencil = np.array( [ - [ - 0.25 * (dx**-2 + dy**-2) * w, - dy**-2 * (1 - w), - 0.25 * (dx**-2 + dy**-2) * w, - ], - [ - dx**-2 * (1 - w), - dx**-2 * dy**-2 * (dx**2 + dy**2) * (w - 2), - dx**-2 * (1 - w), - ], - [ - 0.25 * (dx**-2 + dy**-2) * w, - dy**-2 * (1 - w), - 0.25 * (dx**-2 + dy**-2) * w, - ], + [0.25 * dm2 * w, dxm2 * (1 - w), 0.25 * dm2 * w], + [dym2 * (1 - w), (dxm2 + dym2) * (w - 2), dym2 * (1 - w)], + [0.25 * dm2 * w, dxm2 * (1 - w), 0.25 * dm2 * w], ] ) diff --git a/pde/tools/config.py b/pde/tools/config.py index 3f174f14..a00fc380 100644 --- a/pde/tools/config.py +++ b/pde/tools/config.py @@ -51,8 +51,8 @@ "operators.cartesian.laplacian_2d_corner_weight", 0.0, float, - "Weighting factor for the corner points of the 2d cartesian Laplacian. The " - "standard value is zero, which corresponds to the traditional 5-point stencil. " + "Weighting factor for the corner points of the 2d cartesian Laplacian stencil. " + "The standard value is zero, corresponding to the traditional 5-point stencil. " "Alternative choices are 1/2 (Oono-Puri stencil) and 1/3 (Patra-Karttunen or " "Mehrstellen stencil); see https://en.wikipedia.org/wiki/Nine-point_stencil.", ), diff --git a/tests/grids/test_cartesian_grids.py b/tests/grids/test_cartesian_grids.py index 65f25635..0bc44ad1 100644 --- a/tests/grids/test_cartesian_grids.py +++ b/tests/grids/test_cartesian_grids.py @@ -8,7 +8,7 @@ import numpy as np import pytest -from pde import CartesianGrid, UnitGrid +from pde import CartesianGrid, ScalarField, UnitGrid from pde.grids.boundaries import BoundariesBase, PeriodicityError from pde.grids.operators.common import make_derivative @@ -347,3 +347,15 @@ def test_boundary_coordinates(): np.testing.assert_allclose(c, [[2.0, 0.5], [2.0, 1.5]]) c = grid._boundary_coordinates(axis=0, upper=True, offset=0.5) np.testing.assert_allclose(c, [[1.5, 0.5], [1.5, 1.5]]) + + +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("corner_weight", [1e-8, 1 / 3]) +def test_9point_stencil(periodic, corner_weight): + """Test the 9-point Cartesian stencil.""" + grid = CartesianGrid([[-1, 1], [-1, 1]], [17, 17], periodic=periodic) + field = ScalarField.from_expression(grid, "exp(-x**2 - y**2)") + + reference = field.laplace(bc="auto_periodic_neumann") + test = field.laplace(bc="auto_periodic_neumann", corner_weight=corner_weight) + np.testing.assert_allclose(reference.data, test.data, atol=corner_weight / 3)