From 644edbe58f3fa213a3f0cfb877ee3212607ebfae Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 17:16:11 +0000 Subject: [PATCH 1/9] Add accuracy setting in system. --- src/gui/InfoPanel.cpp | 30 ++++++++++++++++++++++++++ src/gui/InfoPanel.hpp | 3 +++ src/readybase/AbstractRD.cpp | 1 + src/readybase/AbstractRD.hpp | 7 ++++++ src/readybase/FormulaOpenCLImageRD.hpp | 2 ++ 5 files changed, 43 insertions(+) diff --git a/src/gui/InfoPanel.cpp b/src/gui/InfoPanel.cpp index 2c1e5683a..80f8863cb 100644 --- a/src/gui/InfoPanel.cpp +++ b/src/gui/InfoPanel.cpp @@ -58,6 +58,8 @@ const wxString InfoPanel::data_type_label = _("Data type"); const wxString InfoPanel::neighborhood_type_label = _("Neighborhood"); const wxString InfoPanel::neighborhood_range_label = _("Neighborhood range"); const wxString InfoPanel::neighborhood_weight_label = _("Neighborhood weight"); +const wxString InfoPanel::accuracy_label = _("Accuracy"); +const wxString InfoPanel::accuracy_labels[3] = { _("Low"), _("Medium"), _("High") }; // ----------------------------------------------------------------------------- @@ -187,6 +189,11 @@ void InfoPanel::Update(const AbstractRD& system) contents += AppendRow(neighborhood_type_label, neighborhood_type_label, system.GetNeighborhoodType() + "-neighbors", false); } + if (system.HasEditableAccuracyOption()) + { + contents += AppendRow(accuracy_label, accuracy_label, accuracy_labels[static_cast(system.GetAccuracy())], true); + } + contents += AppendRow(block_size_label, block_size_label, wxString::Format(wxT("%d x %d x %d"), system.GetBlockSizeX(),system.GetBlockSizeY(),system.GetBlockSizeZ()), system.HasEditableBlockSize()); @@ -599,6 +606,26 @@ void InfoPanel::ChangeDimensions() // ----------------------------------------------------------------------------- +void InfoPanel::ChangeAccuracy() +{ + const AbstractRD::Accuracy old_val = frame->GetCurrentRDSystem().GetAccuracy(); + + wxArrayString choices; + for (const wxString& label : accuracy_labels) + { + choices.Add(label); + } + wxSingleChoiceDialog dlg(this, _("Accuracy:"), _("Select accuracy:"), + choices); + dlg.SetSelection(static_cast(old_val)); + if (dlg.ShowModal() != wxID_OK) return; + const AbstractRD::Accuracy new_val = static_cast(dlg.GetSelection()); + frame->GetCurrentRDSystem().SetAccuracy(new_val); + Update(frame->GetCurrentRDSystem()); +} + +// ----------------------------------------------------------------------------- + void InfoPanel::ChangeBlockSize() { const AbstractRD& sys = frame->GetCurrentRDSystem(); @@ -680,6 +707,9 @@ void InfoPanel::ChangeInfo(const wxString& label) } else if ( label == dimensions_label ) { ChangeDimensions(); + } else if ( label == accuracy_label ) { + ChangeAccuracy(); + } else if ( label == block_size_label ) { ChangeBlockSize(); diff --git a/src/gui/InfoPanel.hpp b/src/gui/InfoPanel.hpp index 8e183db6d..d09f5d09b 100644 --- a/src/gui/InfoPanel.hpp +++ b/src/gui/InfoPanel.hpp @@ -89,6 +89,8 @@ class InfoPanel : public wxPanel static const wxString neighborhood_type_label; static const wxString neighborhood_range_label; static const wxString neighborhood_weight_label; + static const wxString accuracy_label; + static const wxString accuracy_labels[3]; private: @@ -106,6 +108,7 @@ class InfoPanel : public wxPanel void ChangeFormula(); void ChangeDimensions(); void ChangeBlockSize(); + void ChangeAccuracy(); void ChangeWrapOption(); void ChangeDataType(); diff --git a/src/readybase/AbstractRD.cpp b/src/readybase/AbstractRD.cpp index 4ffebe142..1d7ccc383 100644 --- a/src/readybase/AbstractRD.cpp +++ b/src/readybase/AbstractRD.cpp @@ -31,6 +31,7 @@ using namespace std; AbstractRD::AbstractRD(int data_type) : x_spacing_proportion(0.05) , y_spacing_proportion(0.1) + , accuracy(Accuracy::Medium) { this->timesteps_taken = 0; this->need_reload_formula = true; diff --git a/src/readybase/AbstractRD.hpp b/src/readybase/AbstractRD.hpp index a105beefd..efcca1258 100644 --- a/src/readybase/AbstractRD.hpp +++ b/src/readybase/AbstractRD.hpp @@ -133,6 +133,11 @@ class AbstractRD bool GetWrap() const { return this->wrap; } virtual void SetWrap(bool w) { this->wrap = w; } + virtual bool HasEditableAccuracyOption() const { return false; } + enum class Accuracy { Low, Medium, High }; + Accuracy GetAccuracy() const { return this->accuracy; } + virtual void SetAccuracy(Accuracy acc) { this->accuracy = acc; } + /// Retrieve the current 3D object as a vtkPolyData. virtual void GetAsMesh(vtkPolyData *out,const Properties& render_settings) const =0; @@ -226,6 +231,8 @@ class AbstractRD double x_spacing_proportion; /// spatial separation for rendering multiple chemicals, as a proportion of X double y_spacing_proportion; /// spatial separation for rendering multiple chemicals, as a proportion of Y + Accuracy accuracy; + protected: // functions /// Advance the RD system by n timesteps. diff --git a/src/readybase/FormulaOpenCLImageRD.hpp b/src/readybase/FormulaOpenCLImageRD.hpp index f654639f7..f4de80e20 100644 --- a/src/readybase/FormulaOpenCLImageRD.hpp +++ b/src/readybase/FormulaOpenCLImageRD.hpp @@ -42,6 +42,8 @@ class FormulaOpenCLImageRD : public OpenCLImageRD void SetBlockSizeY(int n) override { this->block_size[1] = n; this->need_reload_formula = true; } void SetBlockSizeZ(int n) override { this->block_size[2] = n; this->need_reload_formula = true; } + bool HasEditableAccuracyOption() const override { return true; } + std::string AssembleKernelSourceFromFormula(const std::string& formula) const override; // we override the parameter access functions because changing the parameters requires rewriting the kernel From bea7d5044fff02acaad2b7937897475e0fcf7a1e Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 17:16:49 +0000 Subject: [PATCH 2/9] Pass accuracy setting through to stencil creation. --- src/readybase/FormulaOpenCLImageRD.cpp | 7 ++++--- src/readybase/FormulaOpenCLImageRD.hpp | 2 ++ src/readybase/stencils.cpp | 16 +++++++++++----- src/readybase/stencils.hpp | 5 ++++- 4 files changed, 21 insertions(+), 9 deletions(-) diff --git a/src/readybase/FormulaOpenCLImageRD.cpp b/src/readybase/FormulaOpenCLImageRD.cpp index 4be044fab..85c691ee5 100644 --- a/src/readybase/FormulaOpenCLImageRD.cpp +++ b/src/readybase/FormulaOpenCLImageRD.cpp @@ -65,12 +65,13 @@ struct InputsNeeded { // ------------------------------------------------------------------------- -InputsNeeded DetectInputsNeeded(const string& formula, int num_chemicals, int dimensionality, const int block_size[3]) +InputsNeeded DetectInputsNeeded(const string& formula, int num_chemicals, int dimensionality, const int block_size[3], + const AbstractRD::Accuracy& accuracy) { InputsNeeded inputs_needed; const vector formula_tokens = tokenize_for_keywords(formula); - const vector known_stencils = GetKnownStencils(dimensionality); + const vector known_stencils = GetKnownStencils(dimensionality, accuracy); for (int i = 0; i < num_chemicals; i++) { const string chem = GetChemicalName(i); @@ -367,7 +368,7 @@ string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(const string& formu } const InputsNeeded inputs_needed = DetectInputsNeeded(formula, this->GetNumberOfChemicals(), - this->GetArenaDimensionality(), this->block_size); + this->GetArenaDimensionality(), this->block_size, this->GetAccuracy()); const string indent = " "; const KernelOptions options(this->wrap, indent, this->data_type, full_data_type_string, this->data_type_suffix, this->block_size); diff --git a/src/readybase/FormulaOpenCLImageRD.hpp b/src/readybase/FormulaOpenCLImageRD.hpp index f4de80e20..1b0cd9572 100644 --- a/src/readybase/FormulaOpenCLImageRD.hpp +++ b/src/readybase/FormulaOpenCLImageRD.hpp @@ -43,6 +43,8 @@ class FormulaOpenCLImageRD : public OpenCLImageRD void SetBlockSizeZ(int n) override { this->block_size[2] = n; this->need_reload_formula = true; } bool HasEditableAccuracyOption() const override { return true; } + void SetAccuracy(Accuracy acc) override { this->accuracy = acc; this->need_reload_formula = true; } + std::string AssembleKernelSourceFromFormula(const std::string& formula) const override; diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index d09ea291e..48fe3235c 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -15,7 +15,6 @@ You should have received a copy of the GNU General Public License along with Ready. If not, see . */ -// Local: #include "stencils.hpp" // Stdlib: @@ -364,7 +363,7 @@ Stencil GetGaussianStencil(int dimensionality) // --------------------------------------------------------------------- -Stencil GetLaplacianStencil(int dimensionality) +Stencil GetLaplacianStencil(int dimensionality, const AbstractRD::Accuracy& accuracy) { switch (dimensionality) { @@ -372,7 +371,14 @@ Stencil GetLaplacianStencil(int dimensionality) return StencilFrom1DArray("laplacian", {1,-2,1}, 1, 2, 0); case 2: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - return StencilFrom2DArray("laplacian", {{1,4,1}, {4,-20,4}, {1,4,1}}, 6, 2, 0, 1); // Known under the name "Mehrstellen" + switch (accuracy) + { + case AbstractRD::Accuracy::Low: + return StencilFrom2DArray("laplacian", { {0,1,0}, {1,-4,1}, {0,1,0} }, 1, 2, 0, 1); // anisotropic + case AbstractRD::Accuracy::Medium: + case AbstractRD::Accuracy::High: + return StencilFrom2DArray("laplacian", { {1,4,1}, {4,-20,4}, {1,4,1} }, 6, 2, 0, 1); // Known under the name "Mehrstellen" + } case 3: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. int c1, c2, c3, c4, divisor; @@ -456,7 +462,7 @@ Stencil GetTriLaplacianStencil(int dimensionality) // --------------------------------------------------------------------- -vector GetKnownStencils(int dimensionality) +vector GetKnownStencils(int dimensionality, const AbstractRD::Accuracy& accuracy) { Stencil XGradient = StencilFrom1DArray("x_gradient", { -1, 0, 1 }, 2, 1, 0); // name, stencil, divisor, dx_power, axes //Stencil XGradient = StencilFrom1DArray("x_gradient", { 1, -8, 0, 8, -1 }, 12, 1, 0); // 4th order version @@ -494,7 +500,7 @@ vector GetKnownStencils(int dimensionality) {-1, 0, 1}, {0, 1, 2} }, 1, 0, 0, 1); Stencil Gaussian = GetGaussianStencil(dimensionality); - Stencil Laplacian = GetLaplacianStencil(dimensionality); + Stencil Laplacian = GetLaplacianStencil(dimensionality, accuracy); Stencil BiLaplacian = GetBiLaplacianStencil(dimensionality); Stencil TriLaplacian = GetTriLaplacianStencil(dimensionality); return { XGradient, YGradient, ZGradient, diff --git a/src/readybase/stencils.hpp b/src/readybase/stencils.hpp index cb4f6c53f..4a05b52eb 100644 --- a/src/readybase/stencils.hpp +++ b/src/readybase/stencils.hpp @@ -15,6 +15,9 @@ You should have received a copy of the GNU General Public License along with Ready. If not, see . */ +// Local: +#include "AbstractRD.hpp" + // Stdlib: #include #include @@ -96,6 +99,6 @@ struct AppliedStencil // --------------------------------------------------------------------- -std::vector GetKnownStencils(int dimensionality); +std::vector GetKnownStencils(int dimensionality, const AbstractRD::Accuracy& accuracy); // --------------------------------------------------------------------- From 64610e39ede2b3752618bee29e6bd37b7d34ebe5 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 19:23:10 +0000 Subject: [PATCH 3/9] Accuracy versions for gradient family and Laplacian family. --- src/readybase/stencils.cpp | 162 ++++++++++++++++++++++++++++--------- 1 file changed, 125 insertions(+), 37 deletions(-) diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index 48fe3235c..e60e428a2 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -18,6 +18,7 @@ #include "stencils.hpp" // Stdlib: +#include #include #include #include @@ -285,6 +286,20 @@ Stencil StencilFrom1DArray(const string& label, int const (&arr)[N], int divisor template Stencil StencilFrom2DArray(const string& label, int const (&arr)[M][N], int divisor, int dx_power, int dim1, int dim2) +{ + array, M> arr2; + for (int j = 0; j < M; j++) + { + for (int i = 0; i < N; i++) + { + arr2[j][i] = arr[j][i]; + } + } + return StencilFrom2DArray(label, arr2, divisor, dx_power, dim1, dim2); +} + +template +Stencil StencilFrom2DArray(const string& label, const array, M>& arr, int divisor, int dx_power, int dim1, int dim2) { if (M % 2 != 1 || N % 2 != 1) { @@ -311,6 +326,23 @@ Stencil StencilFrom2DArray(const string& label, int const (&arr)[M][N], int divi template Stencil StencilFrom3DArray(const string& label, int const (&arr)[L][M][N], int divisor, int dx_power, int dim1, int dim2, int dim3) +{ + array, M>, L> arr2; + for (int k = 0; k < L; k++) + { + for (int j = 0; j < M; j++) + { + for (int i = 0; i < N; i++) + { + arr2[k][j][i] = arr[k][j][i]; + } + } + } + return StencilFrom3DArray(label, arr2, divisor, dx_power, dim1, dim2, dim3); +} + +template +Stencil StencilFrom3DArray(const string& label, const array, M>, L>& arr, int divisor, int dx_power, int dim1, int dim2, int dim3) { if (L % 2 != 1 || M % 2 != 1 || N % 2 != 1) { @@ -339,16 +371,55 @@ Stencil StencilFrom3DArray(const string& label, int const (&arr)[L][M][N], int d // --------------------------------------------------------------------- +using Arr3x3 = array, 3>; +using Arr3x3x3 = array, 3>, 3>; +using Arr5x5 = array, 5>; +using Arr5x5x5 = array, 5>, 5>; + +Arr3x3 RotationallySymmetric3x3(int c1, int c2, int c3) +{ + return { { {c1, c2, c1}, + {c2, c3, c2}, + {c1, c2, c1} } }; +} + +Arr3x3x3 RotationallySymmetric3x3x3(int c1, int c2, int c3, int c4) +{ + return { { { { {c1, c2, c1}, {c2, c3, c2}, {c1, c2, c1} } }, + { { {c2, c3, c2}, {c3, c4, c3}, {c2, c3, c2} } }, + { { {c1, c2, c1}, {c2, c3, c2}, {c1, c2, c1} } } } }; +} + +Arr5x5 RotationallySymmetric5x5(int c1, int c2, int c3, int c4, int c5, int c6) +{ + return { { {c1, c2, c3, c2, c1}, + {c2, c4, c5, c4, c2}, + {c3, c5, c6, c5, c3}, + {c2, c4, c5, c4, c2}, + {c1, c2, c3, c2, c1} } }; +} + +Arr5x5x5 RotationallySymmetric5x5x5(int c1, int c2, int c3, int c4, int c5, int c6, int c7, int c8, int c9, int c10) +{ + return { { { { {c1,c2,c3,c2,c1}, {c2,c4,c5,c4,c2}, {c3,c5,c8, c5,c3}, {c2,c4,c5,c4,c2}, {c1,c2,c3,c2,c1} } }, + { { {c2,c4,c5,c4,c2}, {c4,c6,c7,c6,c4}, {c5,c7,c9, c7,c5}, {c4,c6,c7,c6,c4}, {c2,c4,c5,c4,c2} } }, + { { {c3,c5,c8,c5,c3}, {c5,c7,c9,c7,c5}, {c8,c9,c10,c9,c8}, {c5,c7,c9,c7,c5}, {c3,c5,c8,c5,c3} } }, + { { {c2,c4,c5,c4,c2}, {c4,c6,c7,c6,c4}, {c5,c7,c9, c7,c5}, {c4,c6,c7,c6,c4}, {c2,c4,c5,c4,c2} } }, + { { {c1,c2,c3,c2,c1}, {c2,c4,c5,c4,c2}, {c3,c5,c8, c5,c3}, {c2,c4,c5,c4,c2}, {c1,c2,c3,c2,c1} } } } }; +} + +// --------------------------------------------------------------------- + Stencil GetGaussianStencil(int dimensionality) { switch (dimensionality) { case 1: // from OpenCV - return StencilFrom1DArray("gaussian", {1,4,6,4,1}, 16, 0, 0); // is dx_power correct? + return StencilFrom1DArray("gaussian", {1,4,6,4,1}, 16, 0, 0); case 2: // from https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm - return StencilFrom2DArray("gaussian", {{1,4,7,4,1},{4,16,26,16,4},{7,26,41,26,7},{4,16,26,16,4},{1,4,7,4,1}}, 273, 0, 0, 1); // is dx_power correct? + return StencilFrom2DArray("gaussian", RotationallySymmetric5x5(1,4,7,16,26,41), 273, 0, 0, 1); case 3: // see Scripts/Python/convolve.py return StencilFrom3DArray("gaussian", {{{1,4,6,4,1},{4,16,25,16,4},{7,27,44,27,7},{4,16,25,16,4},{1,4,6,4,1}}, @@ -368,25 +439,28 @@ Stencil GetLaplacianStencil(int dimensionality, const AbstractRD::Accuracy& accu switch (dimensionality) { case 1: - return StencilFrom1DArray("laplacian", {1,-2,1}, 1, 2, 0); + switch (accuracy) + { + case AbstractRD::Accuracy::Low: + case AbstractRD::Accuracy::Medium: // 2nd order version + return StencilFrom1DArray("laplacian", { 1, -2, 1 }, 1, 2, 0); + case AbstractRD::Accuracy::High: // 4th order version + return StencilFrom1DArray("laplacian", { -1, 16, -30, 16, -1 }, 12, 2, 0); + } case 2: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. switch (accuracy) { case AbstractRD::Accuracy::Low: - return StencilFrom2DArray("laplacian", { {0,1,0}, {1,-4,1}, {0,1,0} }, 1, 2, 0, 1); // anisotropic - case AbstractRD::Accuracy::Medium: - case AbstractRD::Accuracy::High: - return StencilFrom2DArray("laplacian", { {1,4,1}, {4,-20,4}, {1,4,1} }, 6, 2, 0, 1); // Known under the name "Mehrstellen" + return StencilFrom2DArray("laplacian", RotationallySymmetric3x3(0, 1, -4), 1, 2, 0, 1); // anisotropic + case AbstractRD::Accuracy::Medium: // 2nd order version + return StencilFrom2DArray("laplacian", RotationallySymmetric3x3(1, 4, -20), 6, 2, 0, 1); // Known under the name "Mehrstellen" + case AbstractRD::Accuracy::High: // 4th order version + return StencilFrom2DArray("laplacian", RotationallySymmetric5x5(0, -2, -1, 16, 52, -252), 60, 2, 0, 1); // 21-point stencil } case 3: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - int c1, c2, c3, c4, divisor; - // 27-point stencil: - c1 = 1; c2 = 3; c3 = 14; c4 = -128; divisor = 30; - return StencilFrom3DArray("laplacian", {{{c1,c2,c1}, {c2,c3,c2}, {c1,c2,c1}}, - {{c2,c3,c2}, {c3,c4,c3}, {c2,c3,c2}}, - {{c1,c2,c1}, {c2,c3,c2}, {c1,c2,c1}}}, divisor, 2, 0, 1, 2); + return StencilFrom3DArray("laplacian", RotationallySymmetric3x3x3(1, 3, 14, -128), 30, 2, 0, 1, 2); // 27-point stencil default: throw runtime_error("Internal error: unsupported dimensionality in GetLaplacianStencil"); } @@ -405,18 +479,10 @@ Stencil GetBiLaplacianStencil(int dimensionality) return StencilFrom1DArray("bilaplacian", {1,-4,6,-4,1}, 1, 4, 0); case 2: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - return StencilFrom2DArray("bilaplacian", {{0,1,1,1,0}, {1,-2,-10,-2,1}, {1,-10,36,-10,1}, {1,-2,-10,-2,1}, {0,1,1,1,0}}, 3, 4, 0, 1); + return StencilFrom2DArray("bilaplacian", RotationallySymmetric5x5(0, 1, 1, -2, -10, 36), 3, 4, 0, 1); case 3: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - int c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, divisor; - // 52-point stencil: - c2 = c3 = c5 = c8 = c9 = 0; - c1 = -1; c4 = 10; c6 = -20; c7 = -36; c10 = 360; divisor = 36; - return StencilFrom3DArray("bilaplacian", {{{c1,c2,c3,c2,c1}, {c2,c4,c5,c4,c2}, {c3,c5,c8, c5,c3}, {c2,c4,c5,c4,c2}, {c1,c2,c3,c2,c1}}, - {{c2,c4,c5,c4,c2}, {c4,c6,c7,c6,c4}, {c5,c7,c9, c7,c5}, {c4,c6,c7,c6,c4}, {c2,c4,c5,c4,c2}}, - {{c3,c5,c8,c5,c3}, {c5,c7,c9,c7,c5}, {c8,c9,c10,c9,c8}, {c5,c7,c9,c7,c5}, {c3,c5,c8,c5,c3}}, - {{c2,c4,c5,c4,c2}, {c4,c6,c7,c6,c4}, {c5,c7,c9, c7,c5}, {c4,c6,c7,c6,c4}, {c2,c4,c5,c4,c2}}, - {{c1,c2,c3,c2,c1}, {c2,c4,c5,c4,c2}, {c3,c5,c8, c5,c3}, {c2,c4,c5,c4,c2}, {c1,c2,c3,c2,c1}}}, divisor, 4, 0, 1, 2); + return StencilFrom3DArray("bilaplacian", RotationallySymmetric5x5x5(-1, 0, 0, 10, 0, -20, -36, 0, 0, 360), 36, 4, 0, 1, 2); // 52-point stencil default: throw runtime_error("Internal error: unsupported dimensionality in GetBiLaplacianStencil"); } @@ -435,9 +501,13 @@ Stencil GetTriLaplacianStencil(int dimensionality) return StencilFrom1DArray("trilaplacian", {1,-6,15,-20,15,-6,1}, 1, 6, 0); case 2: // obtained by convolving a 2D Laplacian stencil with a 2D bi-Laplacian stencil - see Scripts/Python/convolve.py - return StencilFrom2DArray("trilaplacian", { {0,1,5,6,5,1,0}, {1,6,-33,-56,-33,6,1}, {5,-33,6,314,6,-33,5}, - {6,-56,314,-888,314,-56,6}, {5,-33,6,314,6,-33,5}, {1,6,-33,-56,-33,6,1}, - {0,1,5,6,5,1,0} }, 18, 6, 0, 1); + return StencilFrom2DArray("trilaplacian", { {0,1,5,6,5,1,0}, + {1,6,-33,-56,-33,6,1}, + {5,-33,6,314,6,-33,5}, + {6,-56,314,-888,314,-56,6}, + {5,-33,6,314,6,-33,5}, + {1,6,-33,-56,-33,6,1}, + {0,1,5,6,5,1,0} }, 18, 6, 0, 1); case 3: // obtained by convolving a 3D Laplacian stencil with a 3D bi-Laplacian stencil - see Scripts/Python/convolve.py return StencilFrom3DArray("trilaplacian", { {{-1,-3,-1,0,-1,-3,-1},{-3,-4,27,20,27,-4,-3},{-1,27,139,60,139,27,-1}, @@ -464,17 +534,35 @@ Stencil GetTriLaplacianStencil(int dimensionality) vector GetKnownStencils(int dimensionality, const AbstractRD::Accuracy& accuracy) { - Stencil XGradient = StencilFrom1DArray("x_gradient", { -1, 0, 1 }, 2, 1, 0); // name, stencil, divisor, dx_power, axes - //Stencil XGradient = StencilFrom1DArray("x_gradient", { 1, -8, 0, 8, -1 }, 12, 1, 0); // 4th order version - Stencil YGradient = StencilFrom1DArray("y_gradient", { -1, 0, 1 }, 2, 1, 1); - Stencil ZGradient = StencilFrom1DArray("z_gradient", { -1, 0, 1 }, 2, 1, 2); - Stencil XDeriv2 = StencilFrom1DArray("x_deriv2", { 1, -2, 1 }, 1, 2, 0); - Stencil YDeriv2 = StencilFrom1DArray("y_deriv2", { 1, -2, 1 }, 1, 2, 1); - Stencil ZDeriv2 = StencilFrom1DArray("z_deriv2", { 1, -2, 1 }, 1, 2, 2); - Stencil XDeriv3 = StencilFrom1DArray("x_deriv3", { -1, 2, 0, -2, 1 }, 2, 3, 0); - //Stencil XDeriv3 = StencilFrom1DArray("x_deriv3", { 1, -8, 13, 0, -13, 8, -1 }, 8, 3, 0); // 4th order version - Stencil YDeriv3 = StencilFrom1DArray("y_deriv3", { -1, 2, 0, -2, 1 }, 2, 3, 1); - Stencil ZDeriv3 = StencilFrom1DArray("z_deriv3", { -1, 2, 0, -2, 1 }, 2, 3, 2); + Stencil XGradient, YGradient, ZGradient; + Stencil XDeriv2, YDeriv2, ZDeriv2; + Stencil XDeriv3, YDeriv3, ZDeriv3; + switch (accuracy) + { + case AbstractRD::Accuracy::Low: + case AbstractRD::Accuracy::Medium: // 2nd order versions + XGradient = StencilFrom1DArray("x_gradient", { -1, 0, 1 }, 2, 1, 0); + YGradient = StencilFrom1DArray("y_gradient", { -1, 0, 1 }, 2, 1, 1); + ZGradient = StencilFrom1DArray("z_gradient", { -1, 0, 1 }, 2, 1, 2); + XDeriv2 = StencilFrom1DArray("x_deriv2", { 1, -2, 1 }, 1, 2, 0); + YDeriv2 = StencilFrom1DArray("y_deriv2", { 1, -2, 1 }, 1, 2, 1); + ZDeriv2 = StencilFrom1DArray("z_deriv2", { 1, -2, 1 }, 1, 2, 2); + XDeriv3 = StencilFrom1DArray("x_deriv3", { -1, 2, 0, -2, 1 }, 2, 3, 0); + YDeriv3 = StencilFrom1DArray("y_deriv3", { -1, 2, 0, -2, 1 }, 2, 3, 1); + ZDeriv3 = StencilFrom1DArray("z_deriv3", { -1, 2, 0, -2, 1 }, 2, 3, 2); + break; + case AbstractRD::Accuracy::High: // 4th order versions + XGradient = StencilFrom1DArray("x_gradient", { 1, -8, 0, 8, -1 }, 12, 1, 0); + YGradient = StencilFrom1DArray("y_gradient", { 1, -8, 0, 8, -1 }, 12, 1, 1); + ZGradient = StencilFrom1DArray("z_gradient", { 1, -8, 0, 8, -1 }, 12, 1, 2); + XDeriv2 = StencilFrom1DArray("x_deriv2", { -1, 16, -30, 16, -1 }, 12, 2, 0); + YDeriv2 = StencilFrom1DArray("y_deriv2", { -1, 16, -30, 16, -1 }, 12, 2, 1); + ZDeriv2 = StencilFrom1DArray("z_deriv2", { -1, 16, -30, 16, -1 }, 12, 2, 2); + XDeriv3 = StencilFrom1DArray("x_deriv3", { 1, -8, 13, 0, -13, 8, -1 }, 8, 3, 0); + YDeriv3 = StencilFrom1DArray("y_deriv3", { 1, -8, 13, 0, -13, 8, -1 }, 8, 3, 1); + ZDeriv3 = StencilFrom1DArray("z_deriv3", { 1, -8, 13, 0, -13, 8, -1 }, 8, 3, 2); + break; + } Stencil SobelN = StencilFrom2DArray("sobelN", { {1, 2, 1}, {0, 0, 0}, {-1, -2, -1} }, 1, 0, 0, 1); From 14e86c65c81cf99159edf63ac3ed897656b3d505 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 20:43:44 +0000 Subject: [PATCH 4/9] Added accuracy as an optional attribute on the formula node in the XML. --- Help/changes.html | 2 ++ Help/formats.html | 1 + Patterns/advection.vti | 4 +++- src/gui/InfoPanel.cpp | 2 +- src/readybase/FormulaOpenCLImageRD.cpp | 16 ++++++++++++++++ 5 files changed, 23 insertions(+), 2 deletions(-) diff --git a/Help/changes.html b/Help/changes.html index 2da448c48..91fa333d2 100644 --- a/Help/changes.html +++ b/Help/changes.html @@ -21,6 +21,8 @@

Changes

  • Formulas now support a_e, a_nw etc. in a more natural way. Instead of accessing the float4 block that is the neighbor of the current block, these keywords return a float4 that has the four values that are neighbors in the specified way. Now we can do e.g. a = a_ne; in a formula to have the whole pattern shift down and left by one pixel. This is a breaking change: formula rules that accessed cells using a_w, a_e, a_nw, a_se, etc. will no longer behave as they did before. +
  • New pattern setting: accuracy. Controls the size of the stencil. Use "low" to run fast, and "high" to run +accurately.
  • Parameter dx is now reserved for controlling the grid spacing of the Gaussian, Laplacian, bi-Laplacian and tri-Laplacian stencils. Formula rules no longer need to write e.g. laplacian_a / (delta_x * delta_x) and can just write laplacian_a since the value of dx is used in the stencil computation.
  • New render setting: colormap, with 11 maps available: "HSV blend", "spectral", "spectral reversed", "inferno", "inferno reversed", "terrain", "terrain reversed", "orange-purple", "purple-orange", "brown-teal", "teal-brown". The old behavior (HSV interpolation between two colors) is still available - choose "HSV blend". diff --git a/Help/formats.html b/Help/formats.html index ccbff3a7a..ba7e902f3 100644 --- a/Help/formats.html +++ b/Help/formats.html @@ -166,6 +166,7 @@

    <formula>

  • block_size_x (optional) : The x component of the dimensions of the spatial unit processed by each kernel call. Default: 4x1x1
  • block_size_y (optional) : The y component.
  • block_size_z (optional) : The z component. +
  • accuracy (optional) : The stencil accuracy to use. "low", "medium" or "high". Default: "medium".

    Contains:

    An OpenCL kernel snippet, where the chemicals are named a, b, c, etc. diff --git a/Patterns/advection.vti b/Patterns/advection.vti index 124783e82..fdd249861 100644 --- a/Patterns/advection.vti +++ b/Patterns/advection.vti @@ -12,7 +12,9 @@ The rate of change is the negative x gradient. The result is that patterns move to the right without changing shape, as if carried by some bulk flow or <i>advection</i>. - (After a while distortions and instabilities appear because we are using only simple numerical methods.) + After a while distortions and instabilities appear because we are using only simple numerical methods. To reduce + the distortions, set the accuracy to 'high' below to use a 5-point stencil. To delay the onset of numerical instability, + use a smaller timestep or change the data type to 'double'. diff --git a/src/gui/InfoPanel.cpp b/src/gui/InfoPanel.cpp index 80f8863cb..4be2d7c9c 100644 --- a/src/gui/InfoPanel.cpp +++ b/src/gui/InfoPanel.cpp @@ -59,7 +59,7 @@ const wxString InfoPanel::neighborhood_type_label = _("Neighborhood"); const wxString InfoPanel::neighborhood_range_label = _("Neighborhood range"); const wxString InfoPanel::neighborhood_weight_label = _("Neighborhood weight"); const wxString InfoPanel::accuracy_label = _("Accuracy"); -const wxString InfoPanel::accuracy_labels[3] = { _("Low"), _("Medium"), _("High") }; +const wxString InfoPanel::accuracy_labels[3] = { _("low"), _("medium"), _("high") }; // ----------------------------------------------------------------------------- diff --git a/src/readybase/FormulaOpenCLImageRD.cpp b/src/readybase/FormulaOpenCLImageRD.cpp index 85c691ee5..233b9a598 100644 --- a/src/readybase/FormulaOpenCLImageRD.cpp +++ b/src/readybase/FormulaOpenCLImageRD.cpp @@ -414,6 +414,20 @@ void FormulaOpenCLImageRD::InitializeFromXML(vtkXMLDataElement *rd, bool &warn_t // number_of_chemicals: read_required_attribute(xml_formula,"number_of_chemicals",this->n_chemicals); + // accuracy + string accuracy_string; + read_optional_attribute(xml_formula, "accuracy", accuracy_string); + if (accuracy_string.size() > 0) + { + const char* accuracy_labels[3] = { "low", "medium", "high" }; + auto it = find(accuracy_labels, accuracy_labels + 3, accuracy_string); + if (it == accuracy_labels + 3) + { + throw std::runtime_error("unknown accuracy attribute: " + accuracy_string); + } + this->SetAccuracy(static_cast(it - accuracy_labels)); + } + string formula = trim_multiline_string(xml_formula->GetCharacterData()); //this->TestFormula(formula); // will throw on error this->SetFormula(formula); // (won't throw yet) @@ -435,6 +449,8 @@ vtkSmartPointer FormulaOpenCLImageRD::GetAsXML(bool generate_ formula->SetIntAttribute("block_size_x", this->block_size[0]); formula->SetIntAttribute("block_size_y", this->block_size[1]); formula->SetIntAttribute("block_size_z", this->block_size[2]); + const char* accuracy_labels[3] = { "low", "medium", "high" }; + formula->SetAttribute("accuracy", accuracy_labels[static_cast(this->accuracy)]); string f = this->GetFormula(); f = ReplaceAllSubstrings(f, "\n", "\n "); // indent the lines formula->SetCharacterData(f.c_str(), (int)f.length()); From 600548260ff962b1339882df6745242ef771a003 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 20:58:45 +0000 Subject: [PATCH 5/9] Added accuracy option to speed suggestions. --- Help/quickstart.html | 2 ++ src/readybase/FormulaOpenCLImageRD.hpp | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Help/quickstart.html b/Help/quickstart.html index a4dc86444..8e3a83156 100644 --- a/Help/quickstart.html +++ b/Help/quickstart.html @@ -71,6 +71,8 @@

    Quick Start

  • Select the fastest OpenCL device - use Action > Select OpenCL Device...
  • Avoid using the same graphics card for rendering, if possible. For example on my laptop I can tell Windows that Ready should use the integrated graphics chip for rendering (Graphics Settings), which leaves the NVidia card free to run OpenCL. +
  • Change the accuracy to low in the Info Pane, if this does not introduce artifacts into the pattern. View the kernel +to see the effect of this setting.
  • Set the data type to float if possible. Using double is typically slower.
  • Try using local memory, with the setting in the Info Pane. This is a recent feature, let us know if it makes a dramatic difference. diff --git a/src/readybase/FormulaOpenCLImageRD.hpp b/src/readybase/FormulaOpenCLImageRD.hpp index 1b0cd9572..828852afd 100644 --- a/src/readybase/FormulaOpenCLImageRD.hpp +++ b/src/readybase/FormulaOpenCLImageRD.hpp @@ -45,7 +45,6 @@ class FormulaOpenCLImageRD : public OpenCLImageRD bool HasEditableAccuracyOption() const override { return true; } void SetAccuracy(Accuracy acc) override { this->accuracy = acc; this->need_reload_formula = true; } - std::string AssembleKernelSourceFromFormula(const std::string& formula) const override; // we override the parameter access functions because changing the parameters requires rewriting the kernel From 678b5d43d4e94f8dad50fcb7ee55b2e854d64e10 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 21:19:07 +0000 Subject: [PATCH 6/9] Fix: Moved functions to fix gcc/clang error. --- src/readybase/stencils.cpp | 52 +++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index e60e428a2..733c100eb 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -284,20 +284,6 @@ Stencil StencilFrom1DArray(const string& label, int const (&arr)[N], int divisor // --------------------------------------------------------------------- -template -Stencil StencilFrom2DArray(const string& label, int const (&arr)[M][N], int divisor, int dx_power, int dim1, int dim2) -{ - array, M> arr2; - for (int j = 0; j < M; j++) - { - for (int i = 0; i < N; i++) - { - arr2[j][i] = arr[j][i]; - } - } - return StencilFrom2DArray(label, arr2, divisor, dx_power, dim1, dim2); -} - template Stencil StencilFrom2DArray(const string& label, const array, M>& arr, int divisor, int dx_power, int dim1, int dim2) { @@ -322,25 +308,22 @@ Stencil StencilFrom2DArray(const string& label, const array, M>& a return stencil; } -// --------------------------------------------------------------------- - -template -Stencil StencilFrom3DArray(const string& label, int const (&arr)[L][M][N], int divisor, int dx_power, int dim1, int dim2, int dim3) +template +Stencil StencilFrom2DArray(const string& label, int const (&arr)[M][N], int divisor, int dx_power, int dim1, int dim2) { - array, M>, L> arr2; - for (int k = 0; k < L; k++) + array, M> arr2; + for (int j = 0; j < M; j++) { - for (int j = 0; j < M; j++) + for (int i = 0; i < N; i++) { - for (int i = 0; i < N; i++) - { - arr2[k][j][i] = arr[k][j][i]; - } + arr2[j][i] = arr[j][i]; } } - return StencilFrom3DArray(label, arr2, divisor, dx_power, dim1, dim2, dim3); + return StencilFrom2DArray(label, arr2, divisor, dx_power, dim1, dim2); } +// --------------------------------------------------------------------- + template Stencil StencilFrom3DArray(const string& label, const array, M>, L>& arr, int divisor, int dx_power, int dim1, int dim2, int dim3) { @@ -369,6 +352,23 @@ Stencil StencilFrom3DArray(const string& label, const array, return stencil; } +template +Stencil StencilFrom3DArray(const string& label, int const (&arr)[L][M][N], int divisor, int dx_power, int dim1, int dim2, int dim3) +{ + array, M>, L> arr2; + for (int k = 0; k < L; k++) + { + for (int j = 0; j < M; j++) + { + for (int i = 0; i < N; i++) + { + arr2[k][j][i] = arr[k][j][i]; + } + } + } + return StencilFrom3DArray(label, arr2, divisor, dx_power, dim1, dim2, dim3); +} + // --------------------------------------------------------------------- using Arr3x3 = array, 3>; From 1d6286d25cb9d343519d6d796c7ad8a7618e442a Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 21:35:13 +0000 Subject: [PATCH 7/9] Template parameter deduction is failing? --- src/readybase/stencils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index 733c100eb..631d62907 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -419,7 +419,7 @@ Stencil GetGaussianStencil(int dimensionality) return StencilFrom1DArray("gaussian", {1,4,6,4,1}, 16, 0, 0); case 2: // from https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm - return StencilFrom2DArray("gaussian", RotationallySymmetric5x5(1,4,7,16,26,41), 273, 0, 0, 1); + return StencilFrom2DArray<5,5>("gaussian", RotationallySymmetric5x5(1,4,7,16,26,41), 273, 0, 0, 1); case 3: // see Scripts/Python/convolve.py return StencilFrom3DArray("gaussian", {{{1,4,6,4,1},{4,16,25,16,4},{7,27,44,27,7},{4,16,25,16,4},{1,4,6,4,1}}, From 999f6039d02c5dc18acd2441b1cb83c5deba3c35 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 21:59:59 +0000 Subject: [PATCH 8/9] Explicit template parameters to help out some compilers. --- src/readybase/stencils.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index a909bc0ea..67908b6b9 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -360,7 +360,7 @@ Stencil StencilFrom2DArray(const string& label, int const (&arr)[M][N], int divi arr2[j][i] = arr[j][i]; } } - return StencilFrom2DArray(label, arr2, divisor, dx_power, dim1, dim2); + return StencilFrom2DArray(label, arr2, divisor, dx_power, dim1, dim2); } // --------------------------------------------------------------------- @@ -407,7 +407,7 @@ Stencil StencilFrom3DArray(const string& label, int const (&arr)[L][M][N], int d } } } - return StencilFrom3DArray(label, arr2, divisor, dx_power, dim1, dim2, dim3); + return StencilFrom3DArray(label, arr2, divisor, dx_power, dim1, dim2, dim3); } // --------------------------------------------------------------------- @@ -493,15 +493,15 @@ Stencil GetLaplacianStencil(int dimensionality, const AbstractRD::Accuracy& accu switch (accuracy) { case AbstractRD::Accuracy::Low: - return StencilFrom2DArray("laplacian", RotationallySymmetric3x3(0, 1, -4), 1, 2, 0, 1); // anisotropic + return StencilFrom2DArray<3,3>("laplacian", RotationallySymmetric3x3(0, 1, -4), 1, 2, 0, 1); // anisotropic case AbstractRD::Accuracy::Medium: // 2nd order version - return StencilFrom2DArray("laplacian", RotationallySymmetric3x3(1, 4, -20), 6, 2, 0, 1); // Known under the name "Mehrstellen" + return StencilFrom2DArray<3,3>("laplacian", RotationallySymmetric3x3(1, 4, -20), 6, 2, 0, 1); // Known under the name "Mehrstellen" case AbstractRD::Accuracy::High: // 4th order version - return StencilFrom2DArray("laplacian", RotationallySymmetric5x5(0, -2, -1, 16, 52, -252), 60, 2, 0, 1); // 21-point stencil + return StencilFrom2DArray<5,5>("laplacian", RotationallySymmetric5x5(0, -2, -1, 16, 52, -252), 60, 2, 0, 1); // 21-point stencil } case 3: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - return StencilFrom3DArray("laplacian", RotationallySymmetric3x3x3(1, 3, 14, -128), 30, 2, 0, 1, 2); // 27-point stencil + return StencilFrom3DArray<3,3>("laplacian", RotationallySymmetric3x3x3(1, 3, 14, -128), 30, 2, 0, 1, 2); // 27-point stencil default: throw runtime_error("Internal error: unsupported dimensionality in GetLaplacianStencil"); } @@ -520,10 +520,10 @@ Stencil GetBiLaplacianStencil(int dimensionality) return StencilFrom1DArray("bilaplacian", {1,-4,6,-4,1}, 1, 4, 0); case 2: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - return StencilFrom2DArray("bilaplacian", RotationallySymmetric5x5(0, 1, 1, -2, -10, 36), 3, 4, 0, 1); + return StencilFrom2DArray<5,5>("bilaplacian", RotationallySymmetric5x5(0, 1, 1, -2, -10, 36), 3, 4, 0, 1); case 3: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - return StencilFrom3DArray("bilaplacian", RotationallySymmetric5x5x5(-1, 0, 0, 10, 0, -20, -36, 0, 0, 360), 36, 4, 0, 1, 2); // 52-point stencil + return StencilFrom3DArray<5,5,5>("bilaplacian", RotationallySymmetric5x5x5(-1, 0, 0, 10, 0, -20, -36, 0, 0, 360), 36, 4, 0, 1, 2); // 52-point stencil default: throw runtime_error("Internal error: unsupported dimensionality in GetBiLaplacianStencil"); } From 1bfcde308dedfa7b3a2ea0306251e7dfdc136016 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Thu, 21 Jan 2021 22:11:10 +0000 Subject: [PATCH 9/9] Fix: typo --- src/readybase/stencils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index 67908b6b9..4ff391ffb 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -501,7 +501,7 @@ Stencil GetLaplacianStencil(int dimensionality, const AbstractRD::Accuracy& accu } case 3: // Patra, M. & Karttunen, M. (2006) "Stencils with isotropic discretization error for differential operators" Numerical Methods for Partial Differential Equations, 22. - return StencilFrom3DArray<3,3>("laplacian", RotationallySymmetric3x3x3(1, 3, 14, -128), 30, 2, 0, 1, 2); // 27-point stencil + return StencilFrom3DArray<3,3,3>("laplacian", RotationallySymmetric3x3x3(1, 3, 14, -128), 30, 2, 0, 1, 2); // 27-point stencil default: throw runtime_error("Internal error: unsupported dimensionality in GetLaplacianStencil"); }