From 2dbfa074cfca507600c7cdb0998f43f64cc95381 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Mon, 28 Dec 2020 22:00:18 +0000 Subject: [PATCH 01/20] Use swizzling to make keywords a_e etc. be non-block-aligned float4. This is much nicer for the user. --- src/readybase/FormulaOpenCLImageRD.cpp | 73 +++++++++---------- src/readybase/stencils.cpp | 97 +++++++++++++++++--------- src/readybase/stencils.hpp | 5 +- 3 files changed, 106 insertions(+), 69 deletions(-) diff --git a/src/readybase/FormulaOpenCLImageRD.cpp b/src/readybase/FormulaOpenCLImageRD.cpp index cf2bcf137..fea4c28ca 100644 --- a/src/readybase/FormulaOpenCLImageRD.cpp +++ b/src/readybase/FormulaOpenCLImageRD.cpp @@ -60,58 +60,59 @@ struct KeywordOptions { // ------------------------------------------------------------------------- -string GetIndexString(int val, const string& coord, const string& coord_capital, bool wrap) +void AddStencils_Block411(ostringstream& kernel_source, KeywordOptions& options) { - ostringstream oss; - const string index = "index_" + coord; - if (val == 0) - { - oss << index; - } - else if (wrap) - { - oss << "((" << index << showpos << val << " + " << coord_capital << ") & (" << coord_capital << " - 1))"; - } - else + // collect the float4 input blocks we will need + vector input_blocks; + for (const InputPoint& input_point : options.inputs_needed) { - oss << "min(" << coord_capital << "-1, max(0, " << index << showpos << val << "))"; + if (input_point.point.x % 4 == 0) + { + continue; // this block is already aligned + } + if (input_point.point.x >= 0) + { + input_blocks.push_back({ {input_point.point.x + 4 - (input_point.point.x % 4), input_point.point.y, input_point.point.z}, input_point.chem }); + } + else + { + input_blocks.push_back({ {input_point.point.x - 4 - (input_point.point.x % 4), input_point.point.y, input_point.point.z}, input_point.chem }); + } } - return oss.str(); -} - -void AddStencils_Block411(ostringstream& kernel_source, const KeywordOptions& options) -{ - // retrieve the values needed + options.inputs_needed.insert(input_blocks.begin(), input_blocks.end()); + // retrieve the float4 input blocks needed from global memory for (const InputPoint& input_point : options.inputs_needed) { if (input_point.point.x == 0 && input_point.point.y == 0 && input_point.point.z == 0) { continue; // central cell has already been retrieved } - const string index_x = GetIndexString(input_point.point.x, "x", "X", options.wrap); - const string index_y = GetIndexString(input_point.point.y, "y", "Y", options.wrap); - const string index_z = GetIndexString(input_point.point.z, "z", "Z", options.wrap); - kernel_source << options.indent << "const " << options.data_type_string << "4 " << input_point.GetName() - << " = " << input_point.chem << "_in[X * (Y * " << index_z << " + " << index_y << ") + " << index_x << "];\n"; + if (input_point.point.x % 4 == 0) + { + kernel_source << options.indent << "const " << options.data_type_string << "4 " << input_point.GetCode(options.wrap) << ";\n"; + } + } + // compute the non-block-aligned float4's from the block-aligned ones we have retrieved + for (const InputPoint& input_point : options.inputs_needed) + { + if (input_point.point.x % 4 != 0) + { + // swizzle from the retrieved blocks + kernel_source << options.indent << "const " << options.data_type_string << "4 " << input_point.GetName() + << " = (" << options.data_type_string << "4)" << input_point.GetSwizzled() << ";\n"; + } } // compute the stencils needed for (const AppliedStencil& applied_stencil : options.stencils_needed) { kernel_source << options.indent << "const " << options.data_type_string << "4 " << applied_stencil.GetName() - << " = (" << options.data_type_string << "4)(\n" << options.indent << options.indent; - for (int iSlot = 0; iSlot < 4; iSlot++) + << " = ("; + for (int iStencilPoint = 0; iStencilPoint < applied_stencil.stencil.points.size(); iStencilPoint++) { - for (int iStencilPoint = 0; iStencilPoint < applied_stencil.stencil.points.size(); iStencilPoint++) - { - kernel_source << applied_stencil.stencil.points[iStencilPoint].GetCode(iSlot, applied_stencil.chem); - if (iStencilPoint < applied_stencil.stencil.points.size()-1) - { - kernel_source << " + "; - } - } - if (iSlot < 3) + kernel_source << applied_stencil.stencil.points[iStencilPoint].GetCode(applied_stencil.chem); + if (iStencilPoint < applied_stencil.stencil.points.size()-1) { - kernel_source << ",\n" << options.indent << options.indent; + kernel_source << " + "; } } kernel_source << ") / (" << applied_stencil.stencil.GetDivisorCode() << ");\n"; diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index 397141deb..378d59939 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -58,62 +58,99 @@ string Point::GetName() const // --------------------------------------------------------------------- -int CellSlotToBlock(int x) +std::string StencilPoint::GetCode(const string& chem) const { - if (x > 0) + ostringstream oss; + oss << InputPoint{ point, chem }.GetName(); + if (weight != 1) { - return x / 4; + oss << " * " << weight; } - return -int(ceil(-x / 4.0)); + return oss.str(); + // TODO: collect all the stencil points with the same weight and use a single multiply on them } // --------------------------------------------------------------------- -int CellSlotToSlot(int x) +string InputPoint::GetName() const { - while (x < 0) + ostringstream oss; + oss << chem; + const string dir_name = point.GetName(); + if (!dir_name.empty()) { - x += 4; + oss << "_" << dir_name; } - return x % 4; + return oss.str(); } // --------------------------------------------------------------------- -string OffsetToCode(int iSlot, const Point& point, const string& chem) +string InputPoint::GetSwizzled() const { + // assemble a non-block-aligned float4 for the requested point + InputPoint block_left{ point, chem }; + InputPoint block_right{ point, chem }; + if (point.x >= 0) + { + block_left.point.x -= point.x % 4; + block_right.point.x += 4 - (point.x % 4); + } + else + { + block_left.point.x -= 4 + (point.x % 4); + block_right.point.x -= point.x % 4; + } ostringstream oss; - int source_block_x = CellSlotToBlock(iSlot + point.x); - int iSourceSlot = CellSlotToSlot(iSlot + point.x); - InputPoint sourceBlock{ { source_block_x, point.y, point.z }, chem }; - oss << sourceBlock.GetName() << "." << "xyzw"[iSourceSlot]; + oss << "("; + switch (point.x % 4) + { + case 1: + case -3: + oss << block_left.GetName() << ".yzw, " << block_right.GetName() << ".x"; + break; + case 2: + case -2: + oss << block_left.GetName() << ".zw, " << block_right.GetName() << ".xy"; + break; + case 3: + case -1: + oss << block_left.GetName() << ".w, " << block_right.GetName() << ".xyz"; + break; + } + oss << ")"; return oss.str(); } -// --------------------------------------------------------------------- +// ------------------------------------------------------------------------- -std::string StencilPoint::GetCode(int iSlot, const string& chem) const +string InputPoint::GetCode(bool wrap) const { + const string index_x = GetIndexString(point.x / 4, "x", "X", wrap); + const string index_y = GetIndexString(point.y, "y", "Y", wrap); + const string index_z = GetIndexString(point.z, "z", "Z", wrap); ostringstream oss; - oss << OffsetToCode(iSlot, point, chem); - if (weight != 1) - { - oss << " * " << weight; - } + oss << GetName() << " = " << chem << "_in[X * (Y * " << index_z << " + " << index_y << ") + " << index_x << "]"; return oss.str(); - // TODO: collect all the stencil points with the same weight and use a single multiply on them } -// --------------------------------------------------------------------- +// ------------------------------------------------------------------------- -string InputPoint::GetName() const +string InputPoint::GetIndexString(int val, const string& coord, const string& coord_capital, bool wrap) { ostringstream oss; - oss << chem; - const string dir_name = point.GetName(); - if (!dir_name.empty()) + const string index = "index_" + coord; + if (val == 0) { - oss << "_" << dir_name; + oss << index; + } + else if (wrap) + { + oss << "((" << index << showpos << val << " + " << coord_capital << ") & (" << coord_capital << " - 1))"; + } + else + { + oss << "min(" << coord_capital << "-1, max(0, " << index << showpos << val << "))"; } return oss.str(); } @@ -138,11 +175,7 @@ set AppliedStencil::GetInputPoints_Block411() const set input_points; for (const StencilPoint& stencil_point : stencil.points) { - for (int iSlot = 0; iSlot < 4; iSlot++) - { - Point block_point{ CellSlotToBlock(iSlot + stencil_point.point.x), stencil_point.point.y, stencil_point.point.z }; - input_points.insert({ block_point, chem }); - } + input_points.insert({ stencil_point.point, chem }); } return input_points; } diff --git a/src/readybase/stencils.hpp b/src/readybase/stencils.hpp index dedd41f24..866768e3a 100644 --- a/src/readybase/stencils.hpp +++ b/src/readybase/stencils.hpp @@ -48,7 +48,7 @@ struct StencilPoint Point point; int weight; - std::string GetCode(int iSlot, const std::string& chem) const; + std::string GetCode(const std::string& chem) const; }; // --------------------------------------------------------------------- @@ -71,6 +71,9 @@ struct InputPoint std::string chem; std::string GetName() const; + std::string GetCode(bool wrap) const; + std::string GetSwizzled() const; + static std::string GetIndexString(int val, const std::string& coord, const std::string& coord_capital, bool wrap); friend bool operator<(const InputPoint& a, const InputPoint& b) { From 775be1db288e98bdcc366561058bd53afafddd37 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Mon, 28 Dec 2020 22:14:03 +0000 Subject: [PATCH 02/20] Moved code into AppliedStencil::GetCode(). Fix: increased radius of search for chemical access cell keywords. Improved efficiency by only tokenizing formula once. --- src/readybase/FormulaOpenCLImageRD.cpp | 39 ++++++++++---------------- src/readybase/stencils.cpp | 20 +++++++++++-- src/readybase/stencils.hpp | 1 + 3 files changed, 34 insertions(+), 26 deletions(-) diff --git a/src/readybase/FormulaOpenCLImageRD.cpp b/src/readybase/FormulaOpenCLImageRD.cpp index fea4c28ca..45c9db066 100644 --- a/src/readybase/FormulaOpenCLImageRD.cpp +++ b/src/readybase/FormulaOpenCLImageRD.cpp @@ -99,32 +99,21 @@ void AddStencils_Block411(ostringstream& kernel_source, KeywordOptions& options) { // swizzle from the retrieved blocks kernel_source << options.indent << "const " << options.data_type_string << "4 " << input_point.GetName() - << " = (" << options.data_type_string << "4)" << input_point.GetSwizzled() << ";\n"; + << " = (" << options.data_type_string << "4)(" << input_point.GetSwizzled() << ");\n"; } } // compute the stencils needed for (const AppliedStencil& applied_stencil : options.stencils_needed) { - kernel_source << options.indent << "const " << options.data_type_string << "4 " << applied_stencil.GetName() - << " = ("; - for (int iStencilPoint = 0; iStencilPoint < applied_stencil.stencil.points.size(); iStencilPoint++) - { - kernel_source << applied_stencil.stencil.points[iStencilPoint].GetCode(applied_stencil.chem); - if (iStencilPoint < applied_stencil.stencil.points.size()-1) - { - kernel_source << " + "; - } - } - kernel_source << ") / (" << applied_stencil.stencil.GetDivisorCode() << ");\n"; + kernel_source << options.indent << "const " << options.data_type_string << "4 " << applied_stencil.GetCode() << ";\n"; } } // ------------------------------------------------------------------------- -bool UsingKeyword(const string& formula, const string& keyword) +bool UsingKeyword(const vector& formula_tokens, const string& keyword) { - vector tokens = tokenize_for_keywords(formula); - return find(tokens.begin(), tokens.end(), keyword) != tokens.end(); + return find(formula_tokens.begin(), formula_tokens.end(), keyword) != formula_tokens.end(); // TODO: parse properly: ignore comments, not in string, etc. } @@ -174,6 +163,7 @@ std::string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(std::string fo } kernel_source << "\n"; // search the formula for keywords + const vector formula_tokens = tokenize_for_keywords(this->formula); KeywordOptions options{ this->wrap, indent, this->data_type_string, this->data_type_suffix }; const vector known_stencils = GetKnownStencils(this->GetArenaDimensionality()); map gradient_mag_squared; @@ -182,7 +172,7 @@ std::string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(std::string fo const string chem = GetChemicalName(i); // search for keywords that make use of stencils set dependent_stencils; - if (UsingKeyword(formula, "gradient_mag_squared_" + chem)) + if (UsingKeyword(formula_tokens, "gradient_mag_squared_" + chem)) { gradient_mag_squared[chem] = this->GetArenaDimensionality(); switch (this->GetArenaDimensionality()) @@ -200,7 +190,7 @@ std::string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(std::string fo for (const Stencil& stencil : known_stencils) { const string keyword = stencil.label + "_" + chem; - if (UsingKeyword(formula, keyword) || dependent_stencils.find(keyword) != dependent_stencils.end()) + if (UsingKeyword(formula_tokens, keyword) || dependent_stencils.find(keyword) != dependent_stencils.end()) { AppliedStencil applied_stencil{ stencil, chem }; options.stencils_needed.push_back(applied_stencil); @@ -210,15 +200,16 @@ std::string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(std::string fo } } // search for direct access to neighbors, e.g. "a_nw" - for (int x = -1; x <= 1; x++) + const int MAX_RADIUS = 10; // surely if the user wants something this big they should use a kernel? + for (int x = -MAX_RADIUS; x <= MAX_RADIUS; x++) { - for (int y = -2; y <= 2; y++) + for (int y = -MAX_RADIUS; y <= MAX_RADIUS; y++) { - for (int z = -2; z <= 2; z++) + for (int z = -MAX_RADIUS; z <= MAX_RADIUS; z++) { InputPoint input_point{ {x, y, z}, chem }; const string input_point_neighbor_name = input_point.GetName(); - if (UsingKeyword(formula, input_point_neighbor_name)) + if (UsingKeyword(formula_tokens, input_point_neighbor_name)) { options.inputs_needed.insert(input_point); } @@ -241,17 +232,17 @@ std::string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(std::string fo // add the stencils we found AddStencils_Block411(kernel_source, options); // add x_pos, y_pos, z_pos if being used - if (UsingKeyword(formula, "x_pos")) + if (UsingKeyword(formula_tokens, "x_pos")) { kernel_source << indent << "const " << this->data_type_string << "4 x_pos = (4 * index_x + (" << this->data_type_string << "4)(0.0" << this->data_type_suffix << ", 1.0" << this->data_type_suffix << ", 2.0" << this->data_type_suffix << ", 3.0" << this->data_type_suffix << ")) / (4.0" << this->data_type_suffix << " * X);\n"; } - if (UsingKeyword(formula, "y_pos")) + if (UsingKeyword(formula_tokens, "y_pos")) { kernel_source << indent << "const " << this->data_type_string << "4 y_pos = index_y / (" << this->data_type_string << ")(Y); \n"; } - if (UsingKeyword(formula, "z_pos")) + if (UsingKeyword(formula_tokens, "z_pos")) { kernel_source << indent << "const " << this->data_type_string << "4 z_pos = index_z / (" << this->data_type_string << ")(Z);\n"; } diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index 378d59939..66e580f16 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -102,7 +102,6 @@ string InputPoint::GetSwizzled() const block_right.point.x -= point.x % 4; } ostringstream oss; - oss << "("; switch (point.x % 4) { case 1: @@ -118,7 +117,6 @@ string InputPoint::GetSwizzled() const oss << block_left.GetName() << ".w, " << block_right.GetName() << ".xyz"; break; } - oss << ")"; return oss.str(); } @@ -182,6 +180,24 @@ set AppliedStencil::GetInputPoints_Block411() const // --------------------------------------------------------------------- +string AppliedStencil::GetCode() const +{ + ostringstream oss; + oss << GetName() << " = ("; + for (int iStencilPoint = 0; iStencilPoint < stencil.points.size(); iStencilPoint++) + { + oss << stencil.points[iStencilPoint].GetCode(chem); + if (iStencilPoint < stencil.points.size() - 1) + { + oss << " + "; + } + } + oss << ") / (" << stencil.GetDivisorCode() << ")"; + return oss.str(); +} + +// --------------------------------------------------------------------- + template Stencil StencilFrom1DArray(const string& label, int const (&arr)[N], int divisor, int dx_power, int dim1) { diff --git a/src/readybase/stencils.hpp b/src/readybase/stencils.hpp index 866768e3a..7ab1ecdc7 100644 --- a/src/readybase/stencils.hpp +++ b/src/readybase/stencils.hpp @@ -91,6 +91,7 @@ struct AppliedStencil std::string chem; // e.g. "a" std::string GetName() const { return stencil.label + "_" + chem; } + std::string GetCode() const; std::set GetInputPoints_Block411() const; }; From b690fb1fb40eb12f0399f23a5fa828023038a680 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Mon, 28 Dec 2020 22:55:06 +0000 Subject: [PATCH 03/20] Group weights for neater code. Set Gaussian stencil dx_power to zero. --- src/readybase/FormulaOpenCLImageRD.cpp | 2 +- src/readybase/stencils.cpp | 82 ++++++++++++++++++-------- src/readybase/stencils.hpp | 4 +- 3 files changed, 59 insertions(+), 29 deletions(-) diff --git a/src/readybase/FormulaOpenCLImageRD.cpp b/src/readybase/FormulaOpenCLImageRD.cpp index 45c9db066..26ad079cb 100644 --- a/src/readybase/FormulaOpenCLImageRD.cpp +++ b/src/readybase/FormulaOpenCLImageRD.cpp @@ -89,7 +89,7 @@ void AddStencils_Block411(ostringstream& kernel_source, KeywordOptions& options) } if (input_point.point.x % 4 == 0) { - kernel_source << options.indent << "const " << options.data_type_string << "4 " << input_point.GetCode(options.wrap) << ";\n"; + kernel_source << options.indent << "const " << options.data_type_string << "4 " << input_point.GetDirectAccessCode(options.wrap) << ";\n"; } } // compute the non-block-aligned float4's from the block-aligned ones we have retrieved diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index 66e580f16..3d4f81a67 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -21,6 +21,7 @@ // Stdlib: #include #include +#include #include using namespace std; @@ -58,20 +59,6 @@ string Point::GetName() const // --------------------------------------------------------------------- -std::string StencilPoint::GetCode(const string& chem) const -{ - ostringstream oss; - oss << InputPoint{ point, chem }.GetName(); - if (weight != 1) - { - oss << " * " << weight; - } - return oss.str(); - // TODO: collect all the stencil points with the same weight and use a single multiply on them -} - -// --------------------------------------------------------------------- - string InputPoint::GetName() const { ostringstream oss; @@ -122,8 +109,12 @@ string InputPoint::GetSwizzled() const // ------------------------------------------------------------------------- -string InputPoint::GetCode(bool wrap) const +string InputPoint::GetDirectAccessCode(bool wrap) const { + if (point.x % 4 != 0) + { + throw runtime_error("internal error in GetDirectAccessCode: point.x not divisible by 4"); + } const string index_x = GetIndexString(point.x / 4, "x", "X", wrap); const string index_y = GetIndexString(point.y, "y", "Y", wrap); const string index_z = GetIndexString(point.z, "z", "Z", wrap); @@ -158,10 +149,21 @@ string InputPoint::GetIndexString(int val, const string& coord, const string& co string Stencil::GetDivisorCode() const { ostringstream oss; - oss << divisor; - for (int i = 0; i < dx_power; i++) + if (divisor != 1 || dx_power > 0) { - oss << " * dx"; + if (dx_power > 0) + { + oss << "("; + } + oss << divisor; + for (int i = 0; i < dx_power; i++) + { + oss << " * dx"; + } + if (dx_power > 0) + { + oss << ")"; + } } return oss.str(); } @@ -184,15 +186,45 @@ string AppliedStencil::GetCode() const { ostringstream oss; oss << GetName() << " = ("; - for (int iStencilPoint = 0; iStencilPoint < stencil.points.size(); iStencilPoint++) + map> weights; + for (const StencilPoint& stencil_point : stencil.points) + { + weights[stencil_point.weight].push_back(stencil_point.point); + } + bool is_first_weight_list = true; + for(const pair>& weight_list : weights) { - oss << stencil.points[iStencilPoint].GetCode(chem); - if (iStencilPoint < stencil.points.size() - 1) + if (!is_first_weight_list) { oss << " + "; } + is_first_weight_list = false; + const int weight = weight_list.first; + const vector& points = weight_list.second; + if (weight != 1) + { + oss << weight << " * "; + if (points.size() > 1) + { + oss << "("; + } + } + bool is_first_point = true; + for (const Point& point : points) + { + if (!is_first_point) + { + oss << " + "; + } + is_first_point = false; + oss << InputPoint{ point, chem }.GetName(); + } + if (weight != 1 && points.size() > 1) + { + oss << ")"; + } } - oss << ") / (" << stencil.GetDivisorCode() << ")"; + oss << ") / " << stencil.GetDivisorCode(); return oss.str(); } @@ -282,17 +314,17 @@ Stencil GetGaussianStencil(int dimensionality) { case 1: // from OpenCV - return StencilFrom1DArray("gaussian", {1,4,6,4,1}, 16, 1, 0); // is dx_power correct? + return StencilFrom1DArray("gaussian", {1,4,6,4,1}, 16, 0, 0); // is dx_power correct? 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, 1, 0, 1); // is dx_power correct? + 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? 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}}, {{4,16,25,16,4},{16,62,101,62,16},{26,101,164,101,26},{16,62,101,62,16},{4,16,25,16,4}}, {{7,27,44,27,7},{26,101,164,101,26},{41,159,258,159,41},{26,101,164,101,26},{7,27,44,27,7}}, {{4,16,25,16,4},{16,62,101,62,16},{26,101,164,101,26},{16,62,101,62,16},{4,16,25,16,4}}, - {{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}}}, 4390, 1, 0, 1, 2); // is dx_power correct? + {{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}}}, 4390, 0, 0, 1, 2); // is dx_power correct? default: throw runtime_error("Internal error: unsupported dimensionality in GetGaussianStencil"); } diff --git a/src/readybase/stencils.hpp b/src/readybase/stencils.hpp index 7ab1ecdc7..7054b313c 100644 --- a/src/readybase/stencils.hpp +++ b/src/readybase/stencils.hpp @@ -47,8 +47,6 @@ struct StencilPoint { Point point; int weight; - - std::string GetCode(const std::string& chem) const; }; // --------------------------------------------------------------------- @@ -71,7 +69,7 @@ struct InputPoint std::string chem; std::string GetName() const; - std::string GetCode(bool wrap) const; + std::string GetDirectAccessCode(bool wrap) const; std::string GetSwizzled() const; static std::string GetIndexString(int val, const std::string& coord, const std::string& coord_capital, bool wrap); From 229cd7a456a802391612b0ad2c316a068cac40f2 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Mon, 28 Dec 2020 23:15:19 +0000 Subject: [PATCH 04/20] Simplified x_pos code. --- src/readybase/FormulaOpenCLImageRD.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/readybase/FormulaOpenCLImageRD.cpp b/src/readybase/FormulaOpenCLImageRD.cpp index 26ad079cb..1e4de1cef 100644 --- a/src/readybase/FormulaOpenCLImageRD.cpp +++ b/src/readybase/FormulaOpenCLImageRD.cpp @@ -234,9 +234,9 @@ std::string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(std::string fo // add x_pos, y_pos, z_pos if being used if (UsingKeyword(formula_tokens, "x_pos")) { - kernel_source << indent << "const " << this->data_type_string << "4 x_pos = (4 * index_x + (" << this->data_type_string - << "4)(0.0" << this->data_type_suffix << ", 1.0" << this->data_type_suffix << ", 2.0" << this->data_type_suffix - << ", 3.0" << this->data_type_suffix << ")) / (4.0" << this->data_type_suffix << " * X);\n"; + kernel_source << indent << "const " << this->data_type_string << "4 x_pos = (index_x + (" << this->data_type_string + << "4)(0.0" << this->data_type_suffix << ", 0.25" << this->data_type_suffix << ", 0.5" << this->data_type_suffix + << ", 0.75" << this->data_type_suffix << ")) / X;\n"; } if (UsingKeyword(formula_tokens, "y_pos")) { From 52d97ba5f9cf084b66c6bcda7dbd301c473417e7 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 00:22:06 +0000 Subject: [PATCH 05/20] Fixed compass directions, which were a bit muddled. Fixed empty divisor bug. --- .../Experiments/TihaVonGhyczy/sobel_waves.vti | 8 ++-- src/readybase/stencils.cpp | 48 ++++++++++++------- 2 files changed, 34 insertions(+), 22 deletions(-) diff --git a/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti b/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti index 8066552b2..65853bbf1 100644 --- a/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti +++ b/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti @@ -4,17 +4,17 @@ Ocean waves from the simplest of ideas: - da/dt = k1 * laplacian(a) + k2 * sobel_ne(a) + da/dt = k1 * laplacian(a) + k2 * sobel_sw(a) - where sobel_ne is a <a href="https://en.wikipedia.org/wiki/Sobel_operator">Sobel filter convolution</a> - pointing towards the north-east. + where sobel_sw is a <a href="https://en.wikipedia.org/wiki/Sobel_operator">Sobel filter convolution</a> + with positive values coming from the south-west. 0.1 0.1 0.18 - delta_a = k1 * laplacian_a + k2 * sobel_ne_a; + delta_a = k1 * laplacian_a + k2 * sobel_sw_a; diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index 3d4f81a67..a0f593935 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -33,7 +33,7 @@ string Point::GetName() const if (x != 0 || y != 0 || z != 0) { const int order[3] = { 2, 1, 0 }; // output will be e.g. "ne", "usw", "d" - const std::string dirs[3][2] = { {"e", "w"}, {"s", "n"}, {"u", "d"} }; + const std::string dirs[3][2] = { {"e", "w"}, {"n", "s"}, {"u", "d"} }; for (const int i : order) { if (xyz[i] > 0) @@ -151,6 +151,7 @@ string Stencil::GetDivisorCode() const ostringstream oss; if (divisor != 1 || dx_power > 0) { + oss << " / "; if (dx_power > 0) { oss << "("; @@ -185,7 +186,12 @@ set AppliedStencil::GetInputPoints_Block411() const string AppliedStencil::GetCode() const { ostringstream oss; - oss << GetName() << " = ("; + oss << GetName() << " = "; + const string divisor_code = stencil.GetDivisorCode(); + if (!divisor_code.empty()) + { + oss << "("; + } map> weights; for (const StencilPoint& stencil_point : stencil.points) { @@ -224,7 +230,10 @@ string AppliedStencil::GetCode() const oss << ")"; } } - oss << ") / " << stencil.GetDivisorCode(); + if (!divisor_code.empty()) + { + oss << ")" << divisor_code; + } return oss.str(); } @@ -260,16 +269,16 @@ Stencil StencilFrom2DArray(const string& label, int const (&arr)[M][N], int divi throw runtime_error("Internal error: StencilFrom2DArray takes an odd-sized array"); } Stencil stencil{ label, {}, divisor, dx_power }; - for (int j = 0; j < N; j++) + for (int j = 0; j < M; j++) { - for (int i = 0; i < M; i++) + for (int i = 0; i < N; i++) { - if (arr[i][j] != 0) + if (arr[j][i] != 0) { Point point{ 0, 0, 0 }; - point.xyz[dim1] = i - (M - 1) / 2; - point.xyz[dim2] = j - (N - 1) / 2; - stencil.points.push_back({ point, arr[i][j] }); + point.xyz[dim1] = i - (N - 1) / 2; + point.xyz[dim2] = - j + (M - 1) / 2; // rows are in reading order, heading south, which is negative y + stencil.points.push_back({ point, arr[j][i] }); } } } @@ -286,19 +295,19 @@ Stencil StencilFrom3DArray(const string& label, int const (&arr)[L][M][N], int d throw runtime_error("Internal error: StencilFrom3DArray takes an odd-sized array"); } Stencil stencil{ label, {}, divisor, dx_power }; - for (int k = 0; k < N; k++) + for (int k = 0; k < L; k++) { for (int j = 0; j < M; j++) { - for (int i = 0; i < L; i++) + for (int i = 0; i < N; i++) { - if (arr[i][j][k] != 0) + if (arr[k][j][i] != 0) { Point point{ 0, 0, 0 }; - point.xyz[dim1] = i - (L - 1) / 2; - point.xyz[dim2] = j - (M - 1) / 2; - point.xyz[dim3] = k - (N - 1) / 2; - stencil.points.push_back({ point, arr[i][j][k] }); + point.xyz[dim1] = i - (N - 1) / 2; + point.xyz[dim2] = -j + (M - 1) / 2; // rows are in reading order, heading south, which is negative y + point.xyz[dim3] = k - (L - 1) / 2; + stencil.points.push_back({ point, arr[k][j][i] }); } } } @@ -429,12 +438,15 @@ vector GetKnownStencils(int dimensionality) Stencil XGradient = StencilFrom1DArray("x_gradient", { -1, 0, 1 }, 2, 1, 0); Stencil YGradient = StencilFrom1DArray("y_gradient", { -1, 0, 1 }, 2, 1, 1); Stencil ZGradient = StencilFrom1DArray("z_gradient", { -1, 0, 1 }, 2, 1, 2); - Stencil SobelNE = StencilFrom2DArray("sobel_ne", { {2, 1, 0}, {1, 0, -1}, {0, -1, -2} }, 1, 0, 0, 1); + Stencil SobelW = StencilFrom2DArray("sobel_w", { {1, 0, -1}, {2, 0, -2}, {1, 0, -1} }, 1, 0, 0, 1); + Stencil SobelN = StencilFrom2DArray("sobel_n", { {1, 2, 1}, {0, 0, 0}, {-1, -2, -1} }, 1, 0, 0, 1); + Stencil SobelSW = StencilFrom2DArray("sobel_sw", { {0, -1, -2}, {1, 0, -1}, {2, 1, 0} }, 1, 0, 0, 1); Stencil Gaussian = GetGaussianStencil(dimensionality); Stencil Laplacian = GetLaplacianStencil(dimensionality); Stencil BiLaplacian = GetBiLaplacianStencil(dimensionality); Stencil TriLaplacian = GetTriLaplacianStencil(dimensionality); - return { XGradient, YGradient, ZGradient, Gaussian, Laplacian, BiLaplacian, TriLaplacian, SobelNE }; + return { XGradient, YGradient, ZGradient, Gaussian, Laplacian, BiLaplacian, TriLaplacian, + SobelW, SobelN, SobelSW }; } // --------------------------------------------------------------------- From bed8960b7046c4bf1cc053ab5f4f25125ca84641 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 00:48:14 +0000 Subject: [PATCH 06/20] Fix: using non-block-aligned float4's directly would give an error when the block-aligned float4's they need weren't retrieved by something else in the formula. Fix: digits in direct-access keywords were being ignored. --- src/readybase/FormulaOpenCLImageRD.cpp | 7 +++++++ src/readybase/stencils.cpp | 15 +++++++++++++-- src/readybase/stencils.hpp | 1 + src/readybase/utils.cpp | 4 ++-- 4 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/readybase/FormulaOpenCLImageRD.cpp b/src/readybase/FormulaOpenCLImageRD.cpp index 1e4de1cef..9e3af6d91 100644 --- a/src/readybase/FormulaOpenCLImageRD.cpp +++ b/src/readybase/FormulaOpenCLImageRD.cpp @@ -212,6 +212,13 @@ std::string FormulaOpenCLImageRD::AssembleKernelSourceFromFormula(std::string fo if (UsingKeyword(formula_tokens, input_point_neighbor_name)) { options.inputs_needed.insert(input_point); + if (input_point.point.x % 4 != 0) + { + // add the block-aligned float4's we'll need to provide this non-block-aligned float4 + const pair blocks = input_point.GetAlignedBlocks(); + options.inputs_needed.insert(blocks.first); + options.inputs_needed.insert(blocks.second); + } } } } diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index a0f593935..e66e8b247 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -73,9 +73,9 @@ string InputPoint::GetName() const // --------------------------------------------------------------------- -string InputPoint::GetSwizzled() const +pair InputPoint::GetAlignedBlocks() const { - // assemble a non-block-aligned float4 for the requested point + // return the two block-aligned float4's we'll need to assemble this non-block-aligned float4 InputPoint block_left{ point, chem }; InputPoint block_right{ point, chem }; if (point.x >= 0) @@ -88,6 +88,17 @@ string InputPoint::GetSwizzled() const block_left.point.x -= 4 + (point.x % 4); block_right.point.x -= point.x % 4; } + return make_pair(block_left, block_right); +} + +// --------------------------------------------------------------------- + +string InputPoint::GetSwizzled() const +{ + // assemble a non-block-aligned float4 for the requested point + const pair blocks = GetAlignedBlocks(); + const InputPoint& block_left = blocks.first; + const InputPoint& block_right = blocks.second; ostringstream oss; switch (point.x % 4) { diff --git a/src/readybase/stencils.hpp b/src/readybase/stencils.hpp index 7054b313c..3b9fd9536 100644 --- a/src/readybase/stencils.hpp +++ b/src/readybase/stencils.hpp @@ -71,6 +71,7 @@ struct InputPoint std::string GetName() const; std::string GetDirectAccessCode(bool wrap) const; std::string GetSwizzled() const; + std::pair GetAlignedBlocks() const; static std::string GetIndexString(int val, const std::string& coord, const std::string& coord_capital, bool wrap); friend bool operator<(const InputPoint& a, const InputPoint& b) diff --git a/src/readybase/utils.cpp b/src/readybase/utils.cpp index 78ba04fe5..890408823 100644 --- a/src/readybase/utils.cpp +++ b/src/readybase/utils.cpp @@ -189,12 +189,12 @@ string ReplaceAllSubstrings(string subject, const string& search, const string& vector tokenize_for_keywords(const string& formula) { - // customized tokenize for when searching for keywords in formula rules: want case-sensitive, whole words only + // customized tokenize for when searching for keywords in formula rules: whole words only vector tokens; string token; for (char c : formula) { - if ((c >= 'a' && c <= 'z') || c == '_') + if ((c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z') || c == '_' || (c >= '0' && c <= '9')) { token += c; } From e504d1e3518626647b94dc01734022adf0f3155e Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 01:07:14 +0000 Subject: [PATCH 07/20] Added sobelNW, NE, SE for completeness. --- src/readybase/stencils.cpp | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/src/readybase/stencils.cpp b/src/readybase/stencils.cpp index e66e8b247..bda3c50db 100644 --- a/src/readybase/stencils.cpp +++ b/src/readybase/stencils.cpp @@ -449,15 +449,36 @@ vector GetKnownStencils(int dimensionality) Stencil XGradient = StencilFrom1DArray("x_gradient", { -1, 0, 1 }, 2, 1, 0); Stencil YGradient = StencilFrom1DArray("y_gradient", { -1, 0, 1 }, 2, 1, 1); Stencil ZGradient = StencilFrom1DArray("z_gradient", { -1, 0, 1 }, 2, 1, 2); - Stencil SobelW = StencilFrom2DArray("sobel_w", { {1, 0, -1}, {2, 0, -2}, {1, 0, -1} }, 1, 0, 0, 1); - Stencil SobelN = StencilFrom2DArray("sobel_n", { {1, 2, 1}, {0, 0, 0}, {-1, -2, -1} }, 1, 0, 0, 1); - Stencil SobelSW = StencilFrom2DArray("sobel_sw", { {0, -1, -2}, {1, 0, -1}, {2, 1, 0} }, 1, 0, 0, 1); + Stencil SobelN = StencilFrom2DArray("sobelN", { {1, 2, 1}, + {0, 0, 0}, + {-1, -2, -1} }, 1, 0, 0, 1); + Stencil SobelS = StencilFrom2DArray("sobelS", { {-1, -2, -1}, + {0, 0, 0}, + {1, 2, 1} }, 1, 0, 0, 1); + Stencil SobelE = StencilFrom2DArray("sobelE", { {-1, 0, 1}, + {-2, 0, 2}, + {-1, 0, 1} }, 1, 0, 0, 1); + Stencil SobelW = StencilFrom2DArray("sobelW", { {1, 0, -1}, + {2, 0, -2}, + {1, 0, -1} }, 1, 0, 0, 1); + Stencil SobelNW = StencilFrom2DArray("sobelNW", { {2, 1, 0}, + {1, 0, -1}, + {0, -1, -2} }, 1, 0, 0, 1); + Stencil SobelNE = StencilFrom2DArray("sobelNE", { {0, 1, 2}, + {-1, 0, 1}, + {-2, -1, 0} }, 1, 0, 0, 1); + Stencil SobelSW = StencilFrom2DArray("sobelSW", { {0, -1, -2}, + {1, 0, -1}, + {2, 1, 0} }, 1, 0, 0, 1); + Stencil SobelSE = StencilFrom2DArray("sobelSE", { {-2, -1, 0}, + {-1, 0, 1}, + {0, 1, 2} }, 1, 0, 0, 1); Stencil Gaussian = GetGaussianStencil(dimensionality); Stencil Laplacian = GetLaplacianStencil(dimensionality); Stencil BiLaplacian = GetBiLaplacianStencil(dimensionality); Stencil TriLaplacian = GetTriLaplacianStencil(dimensionality); return { XGradient, YGradient, ZGradient, Gaussian, Laplacian, BiLaplacian, TriLaplacian, - SobelW, SobelN, SobelSW }; + SobelN, SobelS, SobelE, SobelW, SobelNW, SobelNE, SobelSW, SobelSE }; } // --------------------------------------------------------------------- From 72a9cb8aa883e10f0a7882606041bd447b697b56 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 01:08:11 +0000 Subject: [PATCH 08/20] Updated sobel_waves to new sobelSW_a naming. --- Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti b/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti index 65853bbf1..24d4ff0aa 100644 --- a/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti +++ b/Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti @@ -4,9 +4,9 @@ Ocean waves from the simplest of ideas: - da/dt = k1 * laplacian(a) + k2 * sobel_sw(a) + da/dt = k1 * laplacian_a + k2 * sobelSW_a - where sobel_sw is a <a href="https://en.wikipedia.org/wiki/Sobel_operator">Sobel filter convolution</a> + where sobelSW_a is a <a href="https://en.wikipedia.org/wiki/Sobel_operator">Sobel filter convolution</a> with positive values coming from the south-west. @@ -14,7 +14,7 @@ 0.1 0.18 - delta_a = k1 * laplacian_a + k2 * sobel_sw_a; + delta_a = k1 * laplacian_a + k2 * sobelSW_a; From 367f5dfe20807ae7d1643320c6b7c66fd60ed5fb Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 01:13:20 +0000 Subject: [PATCH 09/20] Updated surfing-solitons.vti to used Sobel filters and swizzled cell access. --- .../CornusAmmonis/surfing-solitons.vti | 117 +++++++----------- 1 file changed, 45 insertions(+), 72 deletions(-) diff --git a/Patterns/Experiments/CornusAmmonis/surfing-solitons.vti b/Patterns/Experiments/CornusAmmonis/surfing-solitons.vti index f48033cf0..38813b2b6 100644 --- a/Patterns/Experiments/CornusAmmonis/surfing-solitons.vti +++ b/Patterns/Experiments/CornusAmmonis/surfing-solitons.vti @@ -36,41 +36,14 @@ 0.24 - float4 uax = (float4)( -(-(a_nw.w + a_sw.w) + (a_n.y + a_s.y) - 2.0f * a_w.w + 2.0f * a.y ), -(-(a_n.x + a_s.x) + (a_n.z + a_s.z) - 2.0f * a.x + 2.0f * a.z ), -(-(a_n.y + a_s.y) + (a_n.w + a_s.w) - 2.0f * a.y + 2.0f * a.w ), -(-(a_n.z + a_s.z) + (a_ne.x + a_se.x) - 2.0f * a.z + 2.0f * a_e.x)); - -float4 uay = (float4)( -((a_nw.w + a_n.y ) - (a_sw.w + a_s.y) - 2.0f * a_s.x + 2.0f * a_n.x), -((a_n.x + a_n.z ) - (a_s.x + a_s.z) - 2.0f * a_s.y + 2.0f * a_n.y), -((a_n.y + a_n.w ) - (a_s.y + a_s.w) - 2.0f * a_s.z + 2.0f * a_n.z), -((a_n.z + a_ne.x) - (a_s.z + a_se.x) - 2.0f * a_s.w + 2.0f * a_n.w)); - -float4 ubx = (float4)( -(-(b_nw.w + b_sw.w) + (b_n.y + b_s.y) - 2.0f * b_w.w + 2.0f * b.y ), -(-(b_n.x + b_s.x) + (b_n.z + b_s.z) - 2.0f * b.x + 2.0f * b.z ), -(-(b_n.y + b_s.y) + (b_n.w + b_s.w) - 2.0f * b.y + 2.0f * b.w ), -(-(b_n.z + b_s.z) + (b_ne.x + b_se.x) - 2.0f * b.z + 2.0f * b_e.x)); - -float4 uby = (float4)( -((b_nw.w + b_n.y ) - (b_sw.w + b_s.y) - 2.0f * b_s.x + 2.0f * b_n.x), -((b_n.x + b_n.z ) - (b_s.x + b_s.z) - 2.0f * b_s.y + 2.0f * b_n.y), -((b_n.y + b_n.w ) - (b_s.y + b_s.w) - 2.0f * b_s.z + 2.0f * b_n.z), -((b_n.z + b_ne.x) - (b_s.z + b_se.x) - 2.0f * b_s.w + 2.0f * b_n.w)); - -float4 mm_b = (float4)( -(b_nw.w + b_n.x + b_n.y + b_w.w + b.y + b_sw.w + b_s.x + b_s.y ), -(b_n.x + b_n.y + b_n.z + b.x + b.z + b_s.x + b_s.y + b_s.z ), -(b_n.y + b_n.z + b_n.w + b.y + b.w + b_s.y + b_s.z + b_s.w ), -(b_n.z + b_n.w + b_ne.x + b.z + b_e.x + b_s.z + b_s.w + b_se.x)); - -float4 mm_c = (float4)( -(c_nw.w + c_n.x + c_n.y + c_w.w + c.y + c_sw.w + c_s.x + c_s.y ), -(c_n.x + c_n.y + c_n.z + c.x + c.z + c_s.x + c_s.y + c_s.z ), -(c_n.y + c_n.z + c_n.w + c.y + c.w + c_s.y + c_s.z + c_s.w ), -(c_n.z + c_n.w + c_ne.x + c.z + c_e.x + c_s.z + c_s.w + c_se.x)); +float4 uax = sobelE_a; +float4 uay = sobelN_a; + +float4 ubx = sobelE_b; +float4 uby = sobelN_b; + +float4 mm_b = b_nw + b_n + b_ne + b_w + b_e + b_sw + b_s + b_se; +float4 mm_c = c_nw + c_n + c_ne + c_w + c_e + c_sw + c_s + c_se; float4 dp = pow(pow(uax * ubx - uay * uby, 2) + pow(uay * ubx + uax * uby, 2), 0.5f); @@ -79,114 +52,114 @@ d = (mm_c + c) * 0.1f; delta_a = (D_a * (a * a - b * b)) + aScale * laplacian_b + aConst * b - acS * sign(a) * a * a; delta_b = dps * dp + D_b * a * b - bScale * laplacian_a + bConst * a - bcS * sign(b) * b * b; - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + From 9403561993b504afab3850d908ee7b7fb934640f Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 01:54:23 +0000 Subject: [PATCH 10/20] Updated Cornus Ammonis' patterns to use the new keywords. --- .../mutually-catalytic-pinwheels.vti | 190 ++++++---------- ...lly-catalytic-spots-second-order-sobel.vti | 214 ++++++------------ .../pearson-chained-edge-detectors.vti | 77 ++----- .../CornusAmmonis/tip-splitting-web.vti | 119 ++++------ 4 files changed, 205 insertions(+), 395 deletions(-) diff --git a/Patterns/Experiments/CornusAmmonis/mutually-catalytic-pinwheels.vti b/Patterns/Experiments/CornusAmmonis/mutually-catalytic-pinwheels.vti index 5e7bcb97e..e603643a8 100644 --- a/Patterns/Experiments/CornusAmmonis/mutually-catalytic-pinwheels.vti +++ b/Patterns/Experiments/CornusAmmonis/mutually-catalytic-pinwheels.vti @@ -39,255 +39,197 @@ -0.4 - float4 uax = (float4)( -(-(a_nw.w + a_sw.w) + (a_n.y + a_s.y) - 2.0f * a_w.w + 2.0f * a.y ), -(-(a_n.x + a_s.x) + (a_n.z + a_s.z) - 2.0f * a.x + 2.0f * a.z ), -(-(a_n.y + a_s.y) + (a_n.w + a_s.w) - 2.0f * a.y + 2.0f * a.w ), -(-(a_n.z + a_s.z) + (a_ne.x + a_se.x) - 2.0f * a.z + 2.0f * a_e.x)); - - -float4 uay = (float4)( -((a_nw.w + a_n.y ) - (a_sw.w + a_s.y) - 2.0f * a_s.x + 2.0f * a_n.x), -((a_n.x + a_n.z ) - (a_s.x + a_s.z) - 2.0f * a_s.y + 2.0f * a_n.y), -((a_n.y + a_n.w ) - (a_s.y + a_s.w) - 2.0f * a_s.z + 2.0f * a_n.z), -((a_n.z + a_ne.x) - (a_s.z + a_se.x) - 2.0f * a_s.w + 2.0f * a_n.w)); - - -float4 ubx = (float4)( -(-(b_nw.w + b_sw.w) + (b_n.y + b_s.y) - 2.0f * b_w.w + 2.0f * b.y ), -(-(b_n.x + b_s.x) + (b_n.z + b_s.z) - 2.0f * b.x + 2.0f * b.z ), -(-(b_n.y + b_s.y) + (b_n.w + b_s.w) - 2.0f * b.y + 2.0f * b.w ), -(-(b_n.z + b_s.z) + (b_ne.x + b_se.x) - 2.0f * b.z + 2.0f * b_e.x)); - - -float4 uby = (float4)( -((b_nw.w + b_n.y ) - (b_sw.w + b_s.y) - 2.0f * b_s.x + 2.0f * b_n.x), -((b_n.x + b_n.z ) - (b_s.x + b_s.z) - 2.0f * b_s.y + 2.0f * b_n.y), -((b_n.y + b_n.w ) - (b_s.y + b_s.w) - 2.0f * b_s.z + 2.0f * b_n.z), -((b_n.z + b_ne.x) - (b_s.z + b_se.x) - 2.0f * b_s.w + 2.0f * b_n.w)); - - -float4 ucx = (float4)( -(-(c_nw.w + c_sw.w) + (c_n.y + c_s.y) - 2.0f * c_w.w + 2.0f * c.y ), -(-(c_n.x + c_s.x) + (c_n.z + c_s.z) - 2.0f * c.x + 2.0f * c.z ), -(-(c_n.y + c_s.y) + (c_n.w + c_s.w) - 2.0f * c.y + 2.0f * c.w ), -(-(c_n.z + c_s.z) + (c_ne.x + c_se.x) - 2.0f * c.z + 2.0f * c_e.x)); - - -float4 ucy = (float4)( -((c_nw.w + c_n.y ) - (c_sw.w + c_s.y) - 2.0f * c_s.x + 2.0f * c_n.x), -((c_n.x + c_n.z ) - (c_s.x + c_s.z) - 2.0f * c_s.y + 2.0f * c_n.y), -((c_n.y + c_n.w ) - (c_s.y + c_s.w) - 2.0f * c_s.z + 2.0f * c_n.z), -((c_n.z + c_ne.x) - (c_s.z + c_se.x) - 2.0f * c_s.w + 2.0f * c_n.w)); - - - - -float4 udx = (float4)( -(-(d_nw.w + d_sw.w) + (d_n.y + d_s.y) - 2.0f * d_w.w + 2.0f * d.y ), -(-(d_n.x + d_s.x) + (d_n.z + d_s.z) - 2.0f * d.x + 2.0f * d.z ), -(-(d_n.y + d_s.y) + (d_n.w + d_s.w) - 2.0f * d.y + 2.0f * d.w ), -(-(d_n.z + d_s.z) + (d_ne.x + d_se.x) - 2.0f * d.z + 2.0f * d_e.x)); - - -float4 udy = (float4)( -((d_nw.w + d_n.y ) - (d_sw.w + d_s.y) - 2.0f * d_s.x + 2.0f * d_n.x), -((d_n.x + d_n.z ) - (d_s.x + d_s.z) - 2.0f * d_s.y + 2.0f * d_n.y), -((d_n.y + d_n.w ) - (d_s.y + d_s.w) - 2.0f * d_s.z + 2.0f * d_n.z), -((d_n.z + d_ne.x) - (d_s.z + d_se.x) - 2.0f * d_s.w + 2.0f * d_n.w)); - - - - +float4 ubx = sobelE_b; +float4 uby = sobelN_b; +float4 ucx = sobelE_c; +float4 ucy = sobelN_c; +float4 udx = sobelE_d; +float4 udy = sobelN_d; float4 dp1 = ucx * udx + ucy * udy; float4 dp2 = ubx * ucy - uby * ucx; - - - - delta_a = 2.0f * diff_ab * laplacian_a - a*b*b + F1*(1.0f-a); delta_b = diff_ab * laplacian_b + a*b*b - (F1+k1-d*b_feed)*b; delta_c = (a - 0.5f) * dps2 * dp2 + 2.0f * diff_cd * laplacian_c - c*d*d + F2*(1.0f-c); delta_d = dps1 * dp1 + diff_cd * laplacian_d + c*d*d - (F2+k2-step(0.15f,b)*d_feed)*d; e = 2.0f*b + d*2.0f; - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + diff --git a/Patterns/Experiments/CornusAmmonis/mutually-catalytic-spots-second-order-sobel.vti b/Patterns/Experiments/CornusAmmonis/mutually-catalytic-spots-second-order-sobel.vti index 2c7bcea60..acef7325b 100644 --- a/Patterns/Experiments/CornusAmmonis/mutually-catalytic-spots-second-order-sobel.vti +++ b/Patterns/Experiments/CornusAmmonis/mutually-catalytic-spots-second-order-sobel.vti @@ -48,95 +48,25 @@ 0.1 - float4 uax = (float4)( -(-(a_nw.w + a_sw.w) + (a_n.y + a_s.y) - 2.0f * a_w.w + 2.0f * a.y ), -(-(a_n.x + a_s.x) + (a_n.z + a_s.z) - 2.0f * a.x + 2.0f * a.z ), -(-(a_n.y + a_s.y) + (a_n.w + a_s.w) - 2.0f * a.y + 2.0f * a.w ), -(-(a_n.z + a_s.z) + (a_ne.x + a_se.x) - 2.0f * a.z + 2.0f * a_e.x)); - - -float4 uay = (float4)( -((a_nw.w + a_n.y ) - (a_sw.w + a_s.y) - 2.0f * a_s.x + 2.0f * a_n.x), -((a_n.x + a_n.z ) - (a_s.x + a_s.z) - 2.0f * a_s.y + 2.0f * a_n.y), -((a_n.y + a_n.w ) - (a_s.y + a_s.w) - 2.0f * a_s.z + 2.0f * a_n.z), -((a_n.z + a_ne.x) - (a_s.z + a_se.x) - 2.0f * a_s.w + 2.0f * a_n.w)); - - -float4 uax2 = (float4)( -(a_nw.w + a_sw.w + a_n.y + a_s.y) + 2.0f * (a_w.w + a.y ) - 2.0f * (a_n.x + a_s.x + 2.0f * a.x), -(a_n.x + a_s.x + a_n.z + a_s.z) + 2.0f * (a.x + a.z ) - 2.0f * (a_n.y + a_s.y + 2.0f * a.y), -(a_n.y + a_s.y + a_n.w + a_s.w) + 2.0f * (a.y + a.w ) - 2.0f * (a_n.z + a_s.z + 2.0f * a.z), -(a_n.z + a_s.z + a_ne.x + a_se.x) + 2.0f * (a.z + a_e.x) - 2.0f * (a_n.w + a_s.w + 2.0f * a.w)); - - -float4 uay2 = (float4)( -(a_nw.w + a_sw.w + a_n.y + a_s.y) + 2.0f * (a_n.x + a_s.x) - 2.0f * (a_w.w + a.y + 2.0f * a.x), -(a_n.x + a_s.x + a_n.z + a_s.z) + 2.0f * (a_n.y + a_s.y) - 2.0f * (a.x + a.z + 2.0f * a.y), -(a_n.y + a_s.y + a_n.w + a_s.w) + 2.0f * (a_n.z + a_s.z) - 2.0f * (a.y + a.w + 2.0f * a.z), -(a_n.z + a_s.z + a_ne.x + a_se.x) + 2.0f * (a_n.w + a_s.w) - 2.0f * (a.z + a_e.x + 2.0f * a.w)); - - -float4 uaxy = (float4)( -(- a_nw.w + a_sw.w + a_n.y - a_s.y), -(- a_n.x + a_s.x + a_n.z - a_s.z), -(- a_n.y + a_s.y + a_n.w - a_s.w), -(- a_n.z + a_s.z + a_ne.x - a_se.x)); - - +float4 uax = sobelE_a; +float4 uay = sobelN_a; +float4 uax2 = (a_nw + a_sw + a_ne + a_se) + 2.0f * (a_w + a_e) - 2.0f * (a_n + a_s + 2.0f * a); +float4 uay2 = (a_nw + a_sw + a_ne + a_se) + 2.0f * (a_n + a_s) - 2.0f * (a_w + a_e + 2.0f * a); +float4 uaxy = -a_nw + a_sw + a_ne - a_se; float4 uaxsq = pow(uax, 2.0f); float4 uaysq = pow(uay, 2.0f); float4 dsq = 0; if (!(any(isequal(uaxsq, 0)) || any(isequal(uaysq, 0)))) { - dsq = (uaxsq * uax2 + 2 * uax * uay * uaxy + uaysq * uay2) / (uaxsq * uaysq); + dsq = (uaxsq * uax2 + 2 * uax * uay * uaxy + uaysq * uay2) / (uaxsq * uaysq); } - - -float4 ubx = (float4)( -(-(b_nw.w + b_sw.w) + (b_n.y + b_s.y) - 2.0f * b_w.w + 2.0f * b.y ), -(-(b_n.x + b_s.x) + (b_n.z + b_s.z) - 2.0f * b.x + 2.0f * b.z ), -(-(b_n.y + b_s.y) + (b_n.w + b_s.w) - 2.0f * b.y + 2.0f * b.w ), -(-(b_n.z + b_s.z) + (b_ne.x + b_se.x) - 2.0f * b.z + 2.0f * b_e.x)); - - -float4 uby = (float4)( -((b_nw.w + b_n.y ) - (b_sw.w + b_s.y) - 2.0f * b_s.x + 2.0f * b_n.x), -((b_n.x + b_n.z ) - (b_s.x + b_s.z) - 2.0f * b_s.y + 2.0f * b_n.y), -((b_n.y + b_n.w ) - (b_s.y + b_s.w) - 2.0f * b_s.z + 2.0f * b_n.z), -((b_n.z + b_ne.x) - (b_s.z + b_se.x) - 2.0f * b_s.w + 2.0f * b_n.w)); - - -float4 ucx = (float4)( -(-(c_nw.w + c_sw.w) + (c_n.y + c_s.y) - 2.0f * c_w.w + 2.0f * c.y ), -(-(c_n.x + c_s.x) + (c_n.z + c_s.z) - 2.0f * c.x + 2.0f * c.z ), -(-(c_n.y + c_s.y) + (c_n.w + c_s.w) - 2.0f * c.y + 2.0f * c.w ), -(-(c_n.z + c_s.z) + (c_ne.x + c_se.x) - 2.0f * c.z + 2.0f * c_e.x)); - - -float4 ucy = (float4)( -((c_nw.w + c_n.y ) - (c_sw.w + c_s.y) - 2.0f * c_s.x + 2.0f * c_n.x), -((c_n.x + c_n.z ) - (c_s.x + c_s.z) - 2.0f * c_s.y + 2.0f * c_n.y), -((c_n.y + c_n.w ) - (c_s.y + c_s.w) - 2.0f * c_s.z + 2.0f * c_n.z), -((c_n.z + c_ne.x) - (c_s.z + c_se.x) - 2.0f * c_s.w + 2.0f * c_n.w)); - - -float4 udx = (float4)( -(-(d_nw.w + d_sw.w) + (d_n.y + d_s.y) - 2.0f * d_w.w + 2.0f * d.y ), -(-(d_n.x + d_s.x) + (d_n.z + d_s.z) - 2.0f * d.x + 2.0f * d.z ), -(-(d_n.y + d_s.y) + (d_n.w + d_s.w) - 2.0f * d.y + 2.0f * d.w ), -(-(d_n.z + d_s.z) + (d_ne.x + d_se.x) - 2.0f * d.z + 2.0f * d_e.x)); - - -float4 udy = (float4)( -((d_nw.w + d_n.y ) - (d_sw.w + d_s.y) - 2.0f * d_s.x + 2.0f * d_n.x), -((d_n.x + d_n.z ) - (d_s.x + d_s.z) - 2.0f * d_s.y + 2.0f * d_n.y), -((d_n.y + d_n.w ) - (d_s.y + d_s.w) - 2.0f * d_s.z + 2.0f * d_n.z), -((d_n.z + d_ne.x) - (d_s.z + d_se.x) - 2.0f * d_s.w + 2.0f * d_n.w)); - - +float4 ubx = sobelE_b; +float4 uby = sobelN_b; +float4 ucx = sobelE_c; +float4 ucy = sobelN_c; +float4 udx = sobelE_d; +float4 udy = sobelN_d; float4 dp1 = ucx * udx + ucy * udy; float4 dp2 = ubx * ucy - uby * ucx; - - delta_a = 2.0f * diff_ab * laplacian_a - a*b*b + F1*(1.0f-a); delta_b = diff_ab * laplacian_b + a*b*b - (F1+k1-d*b_feed)*b; delta_c = (a - 0.5f) * dps2 * dp2 + 2.0f * diff_cd * laplacian_c - c*d*d + F2*(1.0f-c); @@ -145,183 +75,183 @@ e = 2.0f*b + d*2.0f; f = vs * sign(dsq) * pow(fabs(dsq), vp) + vo; g = (f * 2) + 0.5f; - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + diff --git a/Patterns/Experiments/CornusAmmonis/pearson-chained-edge-detectors.vti b/Patterns/Experiments/CornusAmmonis/pearson-chained-edge-detectors.vti index 7b114c902..c6bc9b516 100644 --- a/Patterns/Experiments/CornusAmmonis/pearson-chained-edge-detectors.vti +++ b/Patterns/Experiments/CornusAmmonis/pearson-chained-edge-detectors.vti @@ -2,9 +2,11 @@ - Chained edge-detector feedback applied to Pearson (1993). Parameter r_s controls edge-detector isotropy. Parameters ds1, ds2, ds3 control scaling of each edge detector layer in the feedback function. If initial conditions result in a stable system, try resetting the initial conditions to see "worm" patterns. + Chained edge-detector feedback applied to Pearson (1993). Parameter r_s controls edge-detector isotropy. + Parameters ds1, ds2, ds3 control scaling of each edge detector layer in the feedback function. If initial + conditions result in a stable system, try resetting the initial conditions to see "worm" patterns. -John E. Pearson (1993) "<a href="http://jaguar.biologie.hu-berlin.de/~wolfram/pages/seminar_theoretische_biologie_2007/literatur/schaber/Pearson1993Science261.pdf">Complex Patterns in a Simple System</a>" Science 261(5118):189-192. + John E. Pearson (1993) "<a href="http://jaguar.biologie.hu-berlin.de/~wolfram/pages/seminar_theoretische_biologie_2007/literatur/schaber/Pearson1993Science261.pdf">Complex Patterns in a Simple System</a>" Science 261(5118):189-192. @@ -44,60 +46,21 @@ John E. Pearson (1993) "<a href="http://jaguar.biologie.hu-berlin.d 3 - float ratio_s = r_s.x; -float vert_s = ratio_s; -float s_s = 2 + vert_s; - -float edge = evs.x / s_s; -float vert = evs.x * vert_s / s_s; - -float4 uax = (float4)( -(edge * (-(a_nw.w + a_sw.w) + (a_n.y + a_s.y)) + vert * (- a_w.w + a.y )), -(edge * (-(a_n.x + a_s.x) + (a_n.z + a_s.z)) + vert * (- a.x + a.z )), -(edge * (-(a_n.y + a_s.y) + (a_n.w + a_s.w)) + vert * (- a.y + a.w )), -(edge * (-(a_n.z + a_s.z) + (a_ne.x + a_se.x)) + vert * (- a.z + a_e.x))); - -float4 uay = (float4)( -(edge * ((a_nw.w + a_n.y ) - (a_sw.w + a_s.y)) + vert * (- a_s.x + a_n.x)), -(edge * ((a_n.x + a_n.z ) - (a_s.x + a_s.z)) + vert * (- a_s.y + a_n.y)), -(edge * ((a_n.y + a_n.w ) - (a_s.y + a_s.w)) + vert * (- a_s.z + a_n.z)), -(edge * ((a_n.z + a_ne.x) - (a_s.z + a_se.x)) + vert * (- a_s.w + a_n.w))); - -float4 ubx = (float4)( -(edge * (-(b_nw.w + b_sw.w) + (b_n.y + b_s.y)) + vert * (- b_w.w + b.y )), -(edge * (-(b_n.x + b_s.x) + (b_n.z + b_s.z)) + vert * (- b.x + b.z )), -(edge * (-(b_n.y + b_s.y) + (b_n.w + b_s.w)) + vert * (- b.y + b.w )), -(edge * (-(b_n.z + b_s.z) + (b_ne.x + b_se.x)) + vert * (- b.z + b_e.x))); - -float4 uby = (float4)( -(edge * ((b_nw.w + b_n.y ) - (b_sw.w + b_s.y)) + vert * (- b_s.x + b_n.x)), -(edge * ((b_n.x + b_n.z ) - (b_s.x + b_s.z)) + vert * (- b_s.y + b_n.y)), -(edge * ((b_n.y + b_n.w ) - (b_s.y + b_s.w)) + vert * (- b_s.z + b_n.z)), -(edge * ((b_n.z + b_ne.x) - (b_s.z + b_se.x)) + vert * (- b_s.w + b_n.w))); - -float4 ucx = (float4)( -(edge * (-(c_nw.w + c_sw.w) + (c_n.y + c_s.y)) + vert * (- c_w.w + c.y )), -(edge * (-(c_n.x + c_s.x) + (c_n.z + c_s.z)) + vert * (- c.x + c.z )), -(edge * (-(c_n.y + c_s.y) + (c_n.w + c_s.w)) + vert * (- c.y + c.w )), -(edge * (-(c_n.z + c_s.z) + (c_ne.x + c_se.x)) + vert * (- c.z + c_e.x))); - -float4 ucy = (float4)( -(edge * ((c_nw.w + c_n.y ) - (c_sw.w + c_s.y)) + vert * (- c_s.x + c_n.x)), -(edge * ((c_n.x + c_n.z ) - (c_s.x + c_s.z)) + vert * (- c_s.y + c_n.y)), -(edge * ((c_n.y + c_n.w ) - (c_s.y + c_s.w)) + vert * (- c_s.z + c_n.z)), -(edge * ((c_n.z + c_ne.x) - (c_s.z + c_se.x)) + vert * (- c_s.w + c_n.w))); - -float4 udx = (float4)( -(edge * (-(d_nw.w + d_sw.w) + (d_n.y + d_s.y)) + vert * (- d_w.w + d.y )), -(edge * (-(d_n.x + d_s.x) + (d_n.z + d_s.z)) + vert * (- d.x + d.z )), -(edge * (-(d_n.y + d_s.y) + (d_n.w + d_s.w)) + vert * (- d.y + d.w )), -(edge * (-(d_n.z + d_s.z) + (d_ne.x + d_se.x)) + vert * (- d.z + d_e.x))); - -float4 udy = (float4)( -(edge * ((d_nw.w + d_n.y ) - (d_sw.w + d_s.y)) + vert * (- d_s.x + d_n.x)), -(edge * ((d_n.x + d_n.z ) - (d_s.x + d_s.z)) + vert * (- d_s.y + d_n.y)), -(edge * ((d_n.y + d_n.w ) - (d_s.y + d_s.w)) + vert * (- d_s.z + d_n.z)), -(edge * ((d_n.z + d_ne.x) - (d_s.z + d_se.x)) + vert * (- d_s.w + d_n.w))); +float4 ratio_s = r_s; +float4 vert_s = ratio_s; +float4 s_s = 2 + vert_s; + +float4 edge = evs / s_s; +float4 vert = evs * vert_s / s_s; + +float4 uax = edge * (-(a_nw + a_sw) + (a_ne + a_se)) + vert * (- a_w + a_e); +float4 uay = edge * ((a_nw + a_ne) - (a_sw + a_se)) + vert * (- a_s + a_n); +float4 ubx = edge * (-(b_nw + b_sw) + (b_ne + b_se)) + vert * (- b_w + b_e); +float4 uby = edge * ((b_nw + b_ne) - (b_sw + b_se)) + vert * (- b_s + b_n); +float4 ucx = edge * (-(c_nw + c_sw) + (c_ne + c_se)) + vert * (- c_w + c_e); +float4 ucy = edge * ((c_nw + c_ne) - (c_sw + c_se)) + vert * (- c_s + c_n); +float4 udx = edge * (-(d_nw + d_sw) + (d_ne + d_se)) + vert * (- d_w + d_e); +float4 udy = edge * ((d_nw + d_ne) - (d_sw + d_se)) + vert * (- d_s + d_n); float4 uam = sqrt(uax * uax + uay * uay); float4 ubm = sqrt(ubx * ubx + uby * uby); @@ -177,7 +140,7 @@ f = ds1 * c + ds2 * d + ds3 * e; - + diff --git a/Patterns/Experiments/CornusAmmonis/tip-splitting-web.vti b/Patterns/Experiments/CornusAmmonis/tip-splitting-web.vti index c489658c7..a3e08cfde 100644 --- a/Patterns/Experiments/CornusAmmonis/tip-splitting-web.vti +++ b/Patterns/Experiments/CornusAmmonis/tip-splitting-web.vti @@ -39,163 +39,138 @@ 0.125 - float4 ua = (float4)( -(a_nw.w + a_n.y + a_sw.w + a_s.y + a_n.x + a_w.w + a.y + a_s.x), -(a_n.x + a_n.z + a_s.x + a_s.z + a_n.y + a.x + a.z + a_s.y), -(a_n.y + a_n.w + a_s.y + a_s.w + a_n.z + a.y + a.w + a_s.z), -(a_n.z + a_ne.x + a_s.z + a_se.x + a_n.w + a.z + a_e.x + a_s.w)); - -float4 uax = (float4)( -(-(a_nw.w + a_sw.w) + (a_n.y + a_s.y) - 2.0 * a_w.w + 2.0 * a.y ), -(-(a_n.x + a_s.x) + (a_n.z + a_s.z) - 2.0 * a.x + 2.0 * a.z ), -(-(a_n.y + a_s.y) + (a_n.w + a_s.w) - 2.0 * a.y + 2.0 * a.w ), -(-(a_n.z + a_s.z) + (a_ne.x + a_se.x) - 2.0 * a.z + 2.0 * a_e.x)); - -float4 uay = (float4)( -((a_nw.w + a_n.y ) - (a_sw.w + a_s.y) - 2.0 * a_s.x + 2.0 * a_n.x), -((a_n.x + a_n.z ) - (a_s.x + a_s.z) - 2.0 * a_s.y + 2.0 * a_n.y), -((a_n.y + a_n.w ) - (a_s.y + a_s.w) - 2.0 * a_s.z + 2.0 * a_n.z), -((a_n.z + a_ne.x) - (a_s.z + a_se.x) - 2.0 * a_s.w + 2.0 * a_n.w)); - -float4 ubx = (float4)( -(-(b_nw.w + b_sw.w) + (b_n.y + b_s.y) - 2.0 * b_w.w + 2.0 * b.y ), -(-(b_n.x + b_s.x) + (b_n.z + b_s.z) - 2.0 * b.x + 2.0 * b.z ), -(-(b_n.y + b_s.y) + (b_n.w + b_s.w) - 2.0 * b.y + 2.0 * b.w ), -(-(b_n.z + b_s.z) + (b_ne.x + b_se.x) - 2.0 * b.z + 2.0 * b_e.x)); - -float4 uby = (float4)( -((b_nw.w + b_n.y ) - (b_sw.w + b_s.y) - 2.0 * b_s.x + 2.0 * b_n.x), -((b_n.x + b_n.z ) - (b_s.x + b_s.z) - 2.0 * b_s.y + 2.0 * b_n.y), -((b_n.y + b_n.w ) - (b_s.y + b_s.w) - 2.0 * b_s.z + 2.0 * b_n.z), -((b_n.z + b_ne.x) - (b_s.z + b_se.x) - 2.0 * b_s.w + 2.0 * b_n.w)); - +float4 ua = a_nw + a_ne + a_sw + a_se + a_n + a_w + a_e + a_s; +float4 uax = sobelE_a; +float4 uay = sobelN_a; +float4 ubx = sobelE_b; +float4 uby = sobelN_b; float4 dp = uax * ubx + uay * uby; delta_a = dpS * dp + k1*a - k2*a*a - a*a*a - b + laplacian_a; delta_b = D_u * (D_us * ua - a) + epsilon*(k3*a - a1*b - a0) + delta*laplacian_b; - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + From 7cd20fd1213b45c643f0b730889254f20f747eea Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 02:35:41 +0000 Subject: [PATCH 11/20] Updated help and changes to the new keywords. --- Help/changes.html | 6 ++++-- Help/writing_new_rules.html | 7 +++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/Help/changes.html b/Help/changes.html index df2161712..9e890c599 100644 --- a/Help/changes.html +++ b/Help/changes.html @@ -13,11 +13,13 @@

Changes

Changes in version ??? (released ???)
  • Formulas now support more stencils: bilaplacian_a, trilaplacian_a, gaussian_a, -sobel_ne_a, x_gradient_a, y_gradient_a, z_gradient_a, gradient_mag_squared_a -(and similarly for the other chemicals). Several patterns have been converted from kernel rules to formula rules +sobelN_a, sobelNE_a, x_gradient_a, gradient_mag_squared_a +(and similarly for the other chemicals and directions). Several patterns have been converted from kernel rules to formula rules as a result.
  • Formulas now support x_pos, y_pos and z_pos, giving the cell location in each direction, in the range [0,1]. This allows a neater and more accurate way of implementing parameter maps. +
  • 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.
  • 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.
  • Moved View Full Kernel... and Show OpenCL Diagnostics... to the View menu. diff --git a/Help/writing_new_rules.html b/Help/writing_new_rules.html index 67625b64a..06fe294f1 100644 --- a/Help/writing_new_rules.html +++ b/Help/writing_new_rules.html @@ -17,7 +17,7 @@

    Formulas and kernel rules

    bilaplacian_aA bi-Laplacian kernel applied to chemical 'a' (or 'b', etc.). Equivalent to applying the Laplacian kernel twice. trilaplacian_aA tri-Laplacian kernel applied to chemical 'a' (or 'b', etc.). Equivalent to applying the Laplacian kernel three times. gaussian_aA Gaussian kernel of sigma=1 applied to chemical 'a' (or 'b', etc.). -sobel_ne_aA Sobel filter applied to chemical 'a' (or 'b', etc.) pointing in the north-east direction. +sobelN_a, sobelE_a, sobelNE_a, etc.A Sobel filter applied to chemical 'a' (or 'b', etc.) pointing in the compass directions. x_gradient_aThe gradient of chemical 'a' (or 'b', etc.) in the x-direction. y_gradient_aThe gradient of chemical 'a' (or 'b', etc.) in the y-direction. z_gradient_aThe gradient of chemical 'a' (or 'b', etc.) in the z-direction. @@ -25,11 +25,14 @@

    Formulas and kernel rules

    x_posThe location of the cell in the x-direction, in the range 0 to 1. y_posThe location of the cell in the y-direction, in the range 0 to 1. z_posThe location of the cell in the z-direction, in the range 0 to 1. +a_n, a_ne, a_n2, a_une, etc.The neighboring cells. Indexed as u=up/d=down Z cells, n=north/s=south Y cells, e=east/w=west X cells: a_[u/d][Z][n/s][Y][e/w][X]. The digit can be omitted if it is one. The digit and the direction can be omitted if the digit is zero.

    The grid spacing (sometimes denoted by h) of the Laplacian and Gaussian stencils is controlled by parameter dx. If there is no parameter called dx then the default value of 1 is used.

    -To understand how the formula is used, try View > View Full Kernel to see the OpenCL kernel code that is assembled. You'll see for example how the timestep parameter appears, at the end during the forward-Euler step. Any formula rule can be permanently converted to a kernel rule by using Action > Convert to Full Kernel. +To understand how the formula is used, try View > View Full Kernel to see the OpenCL kernel code that is assembled. You'll see for example how the timestep parameter is used at the end during the forward Euler step. Any formula rule can be permanently converted to a kernel rule by using Action > Convert to Full Kernel. +

    +For speed we process the image in blocks of 4x1x1, stored as float4. All of the above keywords are float4 values. While doing it this way makes Ready run faster, it does complicate the way we have to write rules - see the examples and please get in touch for more help.

    Images in Ready can be 1D, 2D or 3D. The keywords in the formulas generate stencils of the appropriate dimensionality.

    From 7260d0800825bd425de4caac4d3273ab89a630be Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 02:36:07 +0000 Subject: [PATCH 12/20] Using x_pos in cglrd_ramps_example.vti --- .../DanWills/cglrd_ramps_example.vti | 84 ++++++++++--------- 1 file changed, 44 insertions(+), 40 deletions(-) diff --git a/Patterns/Experiments/DanWills/cglrd_ramps_example.vti b/Patterns/Experiments/DanWills/cglrd_ramps_example.vti index c9dba8b48..63f1724b0 100644 --- a/Patterns/Experiments/DanWills/cglrd_ramps_example.vti +++ b/Patterns/Experiments/DanWills/cglrd_ramps_example.vti @@ -2,11 +2,15 @@ -The complex Ginzburg-Landau equation describes a vast variety of phenomena from nonlinear waves to second-order phase transitions, from superconductivity, superfluidity, and Bose-Einstein condensation to liquid crystals and strings in field theory. +The complex Ginzburg-Landau equation describes a vast variety of phenomena from nonlinear waves to second-order phase +transitions, from superconductivity, superfluidity, and Bose-Einstein condensation to liquid crystals and strings in +field theory. <a href="http://hopf.chem.brandeis.edu/yanglingfa/pattern/cgle/index.html"> [link]</a> -This modification (by Dan Wills) adds a tent-shaped horizontal control ramp to the system, such that the behavior can be made to differ in the center of the frame than at the edges. This is useful as a way to explore the parameter space, as well as allowing the simulation of some interesting ranges of dynamics. +This modification (by Dan Wills) adds a tent-shaped horizontal control ramp to the system, such that the behavior can +be made to differ in the center of the frame than at the edges. This is useful as a way to explore the parameter space, +as well as allowing the simulation of some interesting ranges of dynamics. @@ -52,112 +56,112 @@ This modification (by Dan Wills) adds a tent-shaped horizontal control ramp to t 0 -//make the ramp -d = ( fabs( ( ( get_global_id(0) * 1.0f) / ( get_global_size(0) * 1.0f ) ) - 0.5f ) * 2.0f); -// apply power +// make the ramp +d = fabs( x_pos - 0.5f ) * 2.0f; +// apply power d = pow( d, rampPower ); -//make derivs, applying the ramping +// make derivs, applying the ramping delta_a = (D_a + rampD_a * d) * laplacian_a + (alpha + rampAlpha * d)*a - (gamma + rampGamma * d)*b + (-(beta + rampBeta * d)*a + (delta + rampDelta * d)*b)*(a*a+b*b); delta_b = (D_b + rampD_b * d) * laplacian_b + (alpha + rampAlpha * d)*b + (gamma + rampGamma * d)*a + (-(beta + rampBeta * d)*b - (delta + rampDelta * d)*a)*(a*a+b*b); -//make magnitude +// make magnitude c = sqrt(a*a + b*b); - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + From 512a3a02df61c27710b55d361be4212bec690588 Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 16:56:59 +0000 Subject: [PATCH 13/20] Improved wave_equation.vti --- Patterns/wave_equation.vti | 55 +++++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 22 deletions(-) diff --git a/Patterns/wave_equation.vti b/Patterns/wave_equation.vti index c13169432..083dc9fa6 100644 --- a/Patterns/wave_equation.vti +++ b/Patterns/wave_equation.vti @@ -2,32 +2,43 @@ - Wave equation. Here 'a' stores the value, 'b' stores the rate of change. - - The wave equation says that the second derivative (the rate of change of the rate of change) - is proportional to the Laplacian of the value. Compare with the <a href="open:Patterns/heat_equation.vti">heat equation</a>, where it is - the first derivative (the rate of change) that is proportional to the Laplacian. - - You can reduce the damping parameter to get closer to the pure wave equation but you will need to reduce the timestep too, to avoid - issues with numerical stability. (We use a very simple method: finite differencing with Euler integration.) - - (Parameter 'vis' simply scales the chemicals to have a similar range, for better visualisation.) + The <a href="https://en.wikipedia.org/wiki/Wave_equation">wave equation</a> says that the + second derivative (the rate of change of the rate of change) is proportional to the Laplacian of the value: + + d^2a/dt^2 = laplacian_a + + (Compare with the <a href="open:Patterns/heat_equation.vti">heat equation</a>, where it is + the first derivative (the rate of change) that is proportional to the Laplacian.) + + In our implementation, 'a' stores the value while 'b' stores the rate of change of 'a': + + b = da/dt + + Hence: + + d^2a/dt^2 = db/dt = laplacian_a + + To avoid a time delay in updating 'a' we first update b fully, using our standard forward-Euler integration: + + b += laplacian_a * timestep + + and then we can set the rate of change of 'a': + + delta_a = b - - 0.1 - 0.2 - 15.0 + + 0.6 - delta_a = b/vis + damping*laplacian_a; - delta_b = laplacian_a*vis; + b += laplacian_a * timestep; + delta_a = b; - + - - + + @@ -36,18 +47,18 @@ - + - + - + From 56d15cfb96013b015c518ef45a76d00a393e188a Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 16:57:18 +0000 Subject: [PATCH 14/20] Updated version to 0.11.0 --- CMakeLists.txt | 2 +- Help/about.html | 2 +- Help/changes.html | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 850449311..1a22a17f0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,7 +31,7 @@ endif() project( Ready ) -set( READY_VERSION 0.10.1 ) # check matches Help/about.html +set( READY_VERSION 0.11.0 ) # check matches Help/about.html add_definitions( -D READY_VERSION=${READY_VERSION} ) if( APPLE OR WIN32 ) diff --git a/Help/about.html b/Help/about.html index aef4eaa64..209c3656a 100644 --- a/Help/about.html +++ b/Help/about.html @@ -4,7 +4,7 @@

    -

    Ready, version 0.10.1
    +

    Ready, version 0.11.0
    Copyright 2011-2020 The Ready Bunch:
    diff --git a/Help/changes.html b/Help/changes.html index 9e890c599..10c86747c 100644 --- a/Help/changes.html +++ b/Help/changes.html @@ -10,7 +10,7 @@

    Changes

    https://gollygang.github.io/ready/Help/changes.html.)

    -Changes in version ??? (released ???) +Changes in version 0.11.0 (released ???)

    • Formulas now support more stencils: bilaplacian_a, trilaplacian_a, gaussian_a, sobelN_a, sobelNE_a, x_gradient_a, gradient_mag_squared_a @@ -25,6 +25,7 @@

      Changes

    • Moved View Full Kernel... and Show OpenCL Diagnostics... to the View menu.
    • Mac app requires macOS 10.10 or later.
    • More smaller brush sizes. +
    • Improved wave_equation.vti
    • New patterns:
      • Kobayashi1993/crystals.vti From 67a2625211f0c4120cb3f84c8906a4b0f4e95d6f Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Tue, 29 Dec 2020 17:53:30 +0000 Subject: [PATCH 15/20] Updated smoke.vti to use new keywords for a_n etc. --- Patterns/Experiments/CornusAmmonis/smoke.vti | 237 +++++++------------ 1 file changed, 88 insertions(+), 149 deletions(-) diff --git a/Patterns/Experiments/CornusAmmonis/smoke.vti b/Patterns/Experiments/CornusAmmonis/smoke.vti index 1abf33098..a282f384e 100644 --- a/Patterns/Experiments/CornusAmmonis/smoke.vti +++ b/Patterns/Experiments/CornusAmmonis/smoke.vti @@ -4,251 +4,190 @@ A fast non-physical Eulerian fluid simulation with a wide variety of behaviors. - - - __kernel void rd_compute(__global float *a_in,__global float *b_in,__global float *p_in,__global float *c_in,__global float *a_out,__global float *b_out,__global float *p_out,__global float *c_out) -{ - const float cs = 0.256f; - const float ws = 0.128f; - const float amp = 1.001f; - const float sp = 0.2f; - const float lps = 0.06f; - const float lps2 = -0.08f; - const float ps = 0.0f; - const float sq2 = sqrt(2.0f)/2.0f; - const float wds = 2 + 2 * sqrt(2.0f); - - const int index_x = get_global_id(0); - const int index_y = get_global_id(1); - const int index_z = get_global_id(2); - const int X = get_global_size(0); - const int Y = get_global_size(1); - const int Z = get_global_size(2); - - const int xm1 = ((index_x-1+X) & (X-1)); // wrap (assumes X is a power of 2) - const int xp1 = ((index_x+1) & (X-1)); - const int ym1 = ((index_y-1+Y) & (Y-1)); - const int yp1 = ((index_y+1) & (Y-1)); - const int index_here = X*(Y*index_z + index_y) + index_x; - const int index_n = X*(Y*index_z + ym1) + index_x; - const int index_ne = X*(Y*index_z + ym1) + xp1; - const int index_e = X*(Y*index_z + index_y) + xp1; - const int index_se = X*(Y*index_z + yp1) + xp1; - const int index_s = X*(Y*index_z + yp1) + index_x; - const int index_sw = X*(Y*index_z + yp1) + xm1; - const int index_w = X*(Y*index_z + index_y) + xm1; - const int index_nw = X*(Y*index_z + ym1) + xm1; - - float a = a_in[index_here]; - float b = b_in[index_here]; - float p = p_in[index_here]; - float c = c_in[index_here]; - - float a_n = a_in[index_n]; - float a_ne = a_in[index_ne]; - float a_e = a_in[index_e]; - float a_se = a_in[index_se]; - float a_s = a_in[index_s]; - float a_sw = a_in[index_sw]; - float a_w = a_in[index_w]; - float a_nw = a_in[index_nw]; - - float b_n = b_in[index_n]; - float b_ne = b_in[index_ne]; - float b_e = b_in[index_e]; - float b_se = b_in[index_se]; - float b_s = b_in[index_s]; - float b_sw = b_in[index_sw]; - float b_w = b_in[index_w]; - float b_nw = b_in[index_nw]; - - float p_n = p_in[index_n]; - float p_ne = p_in[index_ne]; - float p_e = p_in[index_e]; - float p_se = p_in[index_se]; - float p_s = p_in[index_s]; - float p_sw = p_in[index_sw]; - float p_w = p_in[index_w]; - float p_nw = p_in[index_nw]; - const float _K0 = -20.0f/6.0f; // center weight - const float _K1 = 4.0f/6.0f; // edge-neighbors - const float _K2 = 1.0f/6.0f; // vertex-neighbors - float laplacian_p = p_n*_K1 + p_ne*_K2 + p_e*_K1 + p_se*_K2 + p_s*_K1 + p_sw*_K2 + p_w*_K1 + p_nw*_K2 + p*_K0; - - float v = (a_n - a_s - b_e + b_w + sq2 * (a_nw + b_nw) + sq2 * (a_ne - b_ne) + sq2 * (b_sw - a_sw) - sq2 * (b_se + a_se)); - float sv = cs * sign(v) * pow(fabs(v), sp); - - float wa = a_w + a_e + sq2 * (a_nw + a_ne + a_se + a_sw) - wds * a; - float wb = b_s + b_n + sq2 * (b_nw + b_ne + b_se + b_sw) - wds * b; - - float pa = a_w + a_e + sq2 * (b_nw - a_nw + a_sw + b_sw + a_ne + b_ne + a_se - b_se); - float pb = b_s + b_n + sq2 * (a_nw - b_nw + a_sw + b_sw + a_ne + b_ne - a_se + b_se); - - p = b_s - b_n - a_e + a_w + sq2 * (a_nw - b_nw) - sq2 * (a_ne + b_ne) + sq2 * (a_sw + b_sw) + sq2 * (b_se - a_se); - - float theta = atan2(b, a); - float mag = sqrt(a*a + b*b); - float cost = cos(theta); - float sint = sin(theta); - float ta = amp * a + ws * wa - lps * cost * laplacian_p + lps2 * a * p; - float tb = amp * b + ws * wb - lps * sint * laplacian_p + lps2 * b * p; - a = ta * cos(sv) - tb * sin(sv); - b = ta * sin(sv) + tb * cos(sv); - - a_out[index_here] = clamp(a, -1.0f, 1.0f); - b_out[index_here] = clamp(b, -1.0f, 1.0f); - p_out[index_here] = p; - - c_out[index_here] = 1000000 * (pb + pa); -} - - + + 0.256 + 0.128 + 1.001 + 0.2 + 0.06 + -0.08 + 0.0 + 0.7071067811865475244 + 4.8284271247461900976 + 1 + +float4 v = (a_n - a_s - b_e + b_w + sq2 * (a_nw + b_nw) + sq2 * (a_ne - b_ne) + sq2 * (b_sw - a_sw) - sq2 * (b_se + a_se)); +float4 sv = cs * sign(v) * pow(fabs(v), sp); + +float4 wa = a_w + a_e + sq2 * (a_nw + a_ne + a_se + a_sw) - wds * a; +float4 wb = b_s + b_n + sq2 * (b_nw + b_ne + b_se + b_sw) - wds * b; + +float4 pa = a_w + a_e + sq2 * (b_nw - a_nw + a_sw + b_sw + a_ne + b_ne + a_se - b_se); +float4 pb = b_s + b_n + sq2 * (a_nw - b_nw + a_sw + b_sw + a_ne + b_ne - a_se + b_se); + +c = b_s - b_n - a_e + a_w + sq2 * (a_nw - b_nw) - sq2 * (a_ne + b_ne) + sq2 * (a_sw + b_sw) + sq2 * (b_se - a_se); + +float4 theta = atan2(b, a); +float4 mag = sqrt(a*a + b*b); +float4 cost = cos(theta); +float4 sint = sin(theta); +float4 ta = amp * a + ws * wa - lps * cost * laplacian_c + lps2 * a * c; +float4 tb = amp * b + ws * wb - lps * sint * laplacian_c + lps2 * b * c; + +a = ta * cos(sv) - tb * sin(sv); +b = ta * sin(sv) + tb * cos(sv); + +a = clamp(a, -1.0f, 1.0f); +b = clamp(b, -1.0f, 1.0f); +d = 1000000 * (pb + pa); + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + From f90a51f4d7e80462b27e93bcecceb58321b2d6c8 Mon Sep 17 00:00:00 2001 From: Dan Wills Date: Wed, 30 Dec 2020 19:46:38 +1030 Subject: [PATCH 16/20] KSE ported to Formula-mode (currently would only work on this branch I believe). Would love it if you could take a look Tim! Pull request owing if it's all-good. --- .../KSE-Formula-stabilize-d3b-256px.vti | 119 ++++++++++++++++++ 1 file changed, 119 insertions(+) create mode 100644 Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti diff --git a/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti b/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti new file mode 100644 index 000000000..db0539efb --- /dev/null +++ b/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti @@ -0,0 +1,119 @@ + + + + + Kuramoto-Sivashinsky-equation, range stabilized (see stabilize parameter) by advice from Steffen Richters. + + Ported to Formula-mode by Dan Wills on the back of Tim Hutton's great work with keyword support. + + It is expected to be black for a while initially (11k timesteps or so) as the values evolve-up to the visualisation-range. + + If you want to reduce dx (to 0.4 say) then you might want to also reduce the timestep (to 0.001, say?) and increase the substeps-per-render to 128 or more, + + + + 0.500000 + + + 0.005000 + + + -0.05 + + + + 0.2 + + + -1.0 + + + -1.0 + + +delta_a = bilapweight * bilaplacian_a + lapweight * laplacian_a + ( gradient_mag_squared_a ) * gms_weight; +delta_a += a * stabilize; +a_out[index_here] = a + timestep * delta_a; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + CAAAAACAAAAAAAAANAAAADQAAAA0AAAANAAAADQAAAA0AAAANAAAADQAAAA=eJztwQEBAAAAgJD+r+4ICgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAYgAAAAXic7cEBAQAAAICQ/q/uCAoAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGIAAAAF4nO3BAQEAAACAkP6v7ggKAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABiAAAABeJztwQEBAAAAgJD+r+4ICgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAYgAAAAXic7cEBAQAAAICQ/q/uCAoAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGIAAAAF4nO3BAQEAAACAkP6v7ggKAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABiAAAABeJztwQEBAAAAgJD+r+4ICgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAYgAAAAXic7cEBAQAAAICQ/q/uCAoAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGIAAAAE= + + + + + + + From 9f3e638bb4aa36599713bb06919bf77272170a72 Mon Sep 17 00:00:00 2001 From: Dan Wills Date: Wed, 30 Dec 2020 20:13:49 +1030 Subject: [PATCH 17/20] Updated Formula-mode VTI for KSE example. --- .../KSE-Formula-stabilize-d3b-256px.vti | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti b/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti index db0539efb..6415b577d 100644 --- a/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti +++ b/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti @@ -2,13 +2,13 @@ - Kuramoto-Sivashinsky-equation, range stabilized (see stabilize parameter) by advice from Steffen Richters. - - Ported to Formula-mode by Dan Wills on the back of Tim Hutton's great work with keyword support. - - It is expected to be black for a while initially (11k timesteps or so) as the values evolve-up to the visualisation-range. - - If you want to reduce dx (to 0.4 say) then you might want to also reduce the timestep (to 0.001, say?) and increase the substeps-per-render to 128 or more, +Kuramoto-Sivashinsky-equation, range stabilized (see the stabilize parameter) based on advice from Steffen Richters ("Richters Finger" on Youtube). + +Ported to Formula-mode by Dan Wills on the back of Tim Hutton's great work with keyword support. + +It is expected to be black for a while initially (11k timesteps or so) as the values evolve-up to the visualisation-range. + +If you want to reduce dx (to 0.4 say) then you might want to also reduce the timestep (to 0.001, say?) and increase the substeps-per-render to 128 or more, @@ -31,9 +31,7 @@ -1.0 -delta_a = bilapweight * bilaplacian_a + lapweight * laplacian_a + ( gradient_mag_squared_a ) * gms_weight; -delta_a += a * stabilize; -a_out[index_here] = a + timestep * delta_a; +delta_a = bilapweight * bilaplacian_a + lapweight * laplacian_a + ( gradient_mag_squared_a ) * gms_weight + a * stabilize; From 5ca1824862af50596c4f946c8521826b35f3496d Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Wed, 30 Dec 2020 10:15:24 +0000 Subject: [PATCH 18/20] Promoted KuramotoSivashinsky1978 and added to CMakeLists.txt and changes.html. Tweaked the initial pattern for faster startup. Tweaked the text. --- CMakeLists.txt | 1 + Help/changes.html | 1 + .../Kuramoto-Sivashinsky_stabilized.vti} | 31 +++++++++---------- 3 files changed, 17 insertions(+), 16 deletions(-) rename Patterns/{Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti => KuramotoSivashinsky1978/Kuramoto-Sivashinsky_stabilized.vti} (80%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1a22a17f0..11e207296 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -262,6 +262,7 @@ set( PATTERN_FILES Patterns/Kobayashi1993/laplacian_growth.vti Patterns/Kobayashi1993/laplacian_growth_3D.vti Patterns/Kobayashi1993/laplacian_growth_3D_corner.vti Patterns/Morozov2008/Fig2.vti Patterns/Morozov2008/Fig4.vti Patterns/Morozov2008/Fig5678_delta_map.vti + Patterns/KuramotoSivashinsky1978/Kuramoto-Sivashinsky_stabilized.vti ) set( HELP_FILES diff --git a/Help/changes.html b/Help/changes.html index 10c86747c..a37a7e670 100644 --- a/Help/changes.html +++ b/Help/changes.html @@ -28,6 +28,7 @@

        Changes

      • Improved wave_equation.vti
      • New patterns:
          +
        • KuramotoSivashinsky1978/Kuramoto-Sivashinsky_stabilized.vti
        • Kobayashi1993/crystals.vti
        • In Experiments/TihaVonGhyczy: sobel_waves.vti.
        • In Experiments/TimHutton: mutually-catalytic_spots_2.vti. diff --git a/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti b/Patterns/KuramotoSivashinsky1978/Kuramoto-Sivashinsky_stabilized.vti similarity index 80% rename from Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti rename to Patterns/KuramotoSivashinsky1978/Kuramoto-Sivashinsky_stabilized.vti index 6415b577d..69ece6fed 100644 --- a/Patterns/Experiments/DanWills/KSE-Formula-stabilize-d3b-256px.vti +++ b/Patterns/KuramotoSivashinsky1978/Kuramoto-Sivashinsky_stabilized.vti @@ -2,43 +2,42 @@ -Kuramoto-Sivashinsky-equation, range stabilized (see the stabilize parameter) based on advice from Steffen Richters ("Richters Finger" on Youtube). +The <a href="https://en.wikipedia.org/wiki/Kuramoto%E2%80%93Sivashinsky_equation">Kuramoto-Sivashinsky-equation</a>, with its range stabilized based on advice from <a href="https://www.youtube.com/channel/UC1Mo9AoX5desRiFLW_skbEg/featured">Steffen Richters</a>. -Ported to Formula-mode by Dan Wills on the back of Tim Hutton's great work with keyword support. +The stabilize parameter prevents the values from growing continuously. You can set it to zero to see what happens without it. -It is expected to be black for a while initially (11k timesteps or so) as the values evolve-up to the visualisation-range. +To increase the size of the blobs, try e.g.: dx = 0.3, timestep = 0.0005 -If you want to reduce dx (to 0.4 say) then you might want to also reduce the timestep (to 0.001, say?) and increase the substeps-per-render to 128 or more, - - +Created by Dan Wills. + + 0.500000 0.005000 - + -0.05 - - - + + 0.2 - + -1.0 - - + + -1.0 - + -delta_a = bilapweight * bilaplacian_a + lapweight * laplacian_a + ( gradient_mag_squared_a ) * gms_weight + a * stabilize; +delta_a = bilapweight * bilaplacian_a + lapweight * laplacian_a + gradient_mag_squared_a * gms_weight + a * stabilize; - + From 2e40ae4328ca87ddee2484f5d3b00621ddf0d47c Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Wed, 30 Dec 2020 11:00:14 +0000 Subject: [PATCH 19/20] Demoted NumericalMethods patterns to Experiments/TimHutton. --- CMakeLists.txt | 40 +++++++++---------- Help/changes.html | 20 +++++----- .../advection_2stepAdamsBashforth.vti | 0 .../advection_2stepAdamsBashforth_double.vti | 0 .../advection_3stepAdamsBashforth.vti | 0 .../advection_3stepAdamsBashforth_double.vti | 0 .../advection_HeunsMethod.vti | 0 .../advection_HeunsMethod_double.vti | 0 .../advection_LaxFriedrichs.vti | 0 .../advection_LaxFriedrichs_double.vti | 0 .../advection_LaxWendroff.vti | 0 .../advection_LaxWendroff_double.vti | 0 .../advection_RungeKutta4.vti | 0 .../advection_RungeKutta4_double.vti | 0 .../advection_forwardEuler.vti | 0 .../advection_forwardEuler_double.vti | 0 .../advection_midpointMethod.vti | 0 .../advection_midpointMethod_double.vti | 0 .../advection_modifiedEuler.vti | 0 .../advection_modifiedEuler_double.vti | 0 .../advection_staggeredLeapfrog.vti | 0 .../advection_staggeredLeapfrog_double.vti | 0 22 files changed, 30 insertions(+), 30 deletions(-) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_2stepAdamsBashforth.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_2stepAdamsBashforth_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_3stepAdamsBashforth.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_3stepAdamsBashforth_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_HeunsMethod.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_HeunsMethod_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_LaxFriedrichs.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_LaxFriedrichs_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_LaxWendroff.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_LaxWendroff_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_RungeKutta4.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_RungeKutta4_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_forwardEuler.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_forwardEuler_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_midpointMethod.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_midpointMethod_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_modifiedEuler.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_modifiedEuler_double.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_staggeredLeapfrog.vti (100%) rename Patterns/{ => Experiments/TimHutton}/NumericalMethods/advection_staggeredLeapfrog_double.vti (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 11e207296..6cf464a94 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -191,6 +191,26 @@ set( PATTERN_FILES Patterns/Experiments/TimHutton/LifeBlur.vti Patterns/Experiments/TimHutton/mutually-catalytic_spots.vti Patterns/Experiments/TimHutton/mutually-catalytic_spots_2.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_forwardEuler.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_2stepAdamsBashforth.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_3stepAdamsBashforth.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_modifiedEuler.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_midpointMethod.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_HeunsMethod.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxFriedrichs.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxWendroff.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_RungeKutta4.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_staggeredLeapfrog.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_forwardEuler_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_2stepAdamsBashforth_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_3stepAdamsBashforth_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_modifiedEuler_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_midpointMethod_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_HeunsMethod_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxFriedrichs_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxWendroff_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_RungeKutta4_double.vti + Patterns/Experiments/TimHutton/NumericalMethods/advection_staggeredLeapfrog_double.vti Patterns/Experiments/DanWills/cglrd_ramps_example.vti Patterns/Experiments/DanWills/grayscott-historyWave_fuseWorms.vti Patterns/Experiments/DanWills/grayscott-historyWave_moreLifelike.vti @@ -230,26 +250,6 @@ set( PATTERN_FILES Patterns/Experiments/TihaVonGhyczy/sobel_waves.vti Patterns/Agmon2014/cells.vti Patterns/Agmon2014/cells_double.vti Patterns/Agmon2014/oil-water.vti - Patterns/NumericalMethods/advection_forwardEuler.vti - Patterns/NumericalMethods/advection_2stepAdamsBashforth.vti - Patterns/NumericalMethods/advection_3stepAdamsBashforth.vti - Patterns/NumericalMethods/advection_modifiedEuler.vti - Patterns/NumericalMethods/advection_midpointMethod.vti - Patterns/NumericalMethods/advection_HeunsMethod.vti - Patterns/NumericalMethods/advection_LaxFriedrichs.vti - Patterns/NumericalMethods/advection_LaxWendroff.vti - Patterns/NumericalMethods/advection_RungeKutta4.vti - Patterns/NumericalMethods/advection_staggeredLeapfrog.vti - Patterns/NumericalMethods/advection_forwardEuler_double.vti - Patterns/NumericalMethods/advection_2stepAdamsBashforth_double.vti - Patterns/NumericalMethods/advection_3stepAdamsBashforth_double.vti - Patterns/NumericalMethods/advection_modifiedEuler_double.vti - Patterns/NumericalMethods/advection_midpointMethod_double.vti - Patterns/NumericalMethods/advection_HeunsMethod_double.vti - Patterns/NumericalMethods/advection_LaxFriedrichs_double.vti - Patterns/NumericalMethods/advection_LaxWendroff_double.vti - Patterns/NumericalMethods/advection_RungeKutta4_double.vti - Patterns/NumericalMethods/advection_staggeredLeapfrog_double.vti Patterns/Guo2014/guo.vti Patterns/Maginu1975/maginu_parameter_map.vti Patterns/Kryuchkov2020/Drosophila_corneal_nanocoatings.vti diff --git a/Help/changes.html b/Help/changes.html index a37a7e670..622473d04 100644 --- a/Help/changes.html +++ b/Help/changes.html @@ -130,16 +130,16 @@

          Changes

        • And more amazing rules by Cornus Ammonis: bz-warpsharp.vti applying the warp-sharp algorithm to the Belousov-Zhabotinsky system, smoke.vti, smoke-ising.vti, and splats.vti.
        • advection.vti: the advection equation, plus some different methods of integrating it, for comparison: - forward Euler, - two-step Adams-Bashforth, - three-step Adams-Bashforth, - modified Euler, - the midpoint method, - Heun's method, - staggered leapfrog, - Lax-Friedrichs, - Lax-Wendroff and - fourth-order Runge-Kutta. + forward Euler, + two-step Adams-Bashforth, + three-step Adams-Bashforth, + modified Euler, + the midpoint method, + Heun's method, + staggered leapfrog, + Lax-Friedrichs, + Lax-Wendroff and + fourth-order Runge-Kutta.
        • Pennybacker2013/phyllotaxis_fibonacci.vti and Pennybacker2013/phyllotaxis_hexagons.vti showing a single-chemical model of phyllotaxis in plants. Also Pennybacker2013/spots.vti, Pennybacker2013/stripes.vti and Pennybacker2013/parameter_map.vti, showing a simpler model without the gradient term. diff --git a/Patterns/NumericalMethods/advection_2stepAdamsBashforth.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_2stepAdamsBashforth.vti similarity index 100% rename from Patterns/NumericalMethods/advection_2stepAdamsBashforth.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_2stepAdamsBashforth.vti diff --git a/Patterns/NumericalMethods/advection_2stepAdamsBashforth_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_2stepAdamsBashforth_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_2stepAdamsBashforth_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_2stepAdamsBashforth_double.vti diff --git a/Patterns/NumericalMethods/advection_3stepAdamsBashforth.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_3stepAdamsBashforth.vti similarity index 100% rename from Patterns/NumericalMethods/advection_3stepAdamsBashforth.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_3stepAdamsBashforth.vti diff --git a/Patterns/NumericalMethods/advection_3stepAdamsBashforth_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_3stepAdamsBashforth_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_3stepAdamsBashforth_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_3stepAdamsBashforth_double.vti diff --git a/Patterns/NumericalMethods/advection_HeunsMethod.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_HeunsMethod.vti similarity index 100% rename from Patterns/NumericalMethods/advection_HeunsMethod.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_HeunsMethod.vti diff --git a/Patterns/NumericalMethods/advection_HeunsMethod_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_HeunsMethod_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_HeunsMethod_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_HeunsMethod_double.vti diff --git a/Patterns/NumericalMethods/advection_LaxFriedrichs.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxFriedrichs.vti similarity index 100% rename from Patterns/NumericalMethods/advection_LaxFriedrichs.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxFriedrichs.vti diff --git a/Patterns/NumericalMethods/advection_LaxFriedrichs_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxFriedrichs_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_LaxFriedrichs_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxFriedrichs_double.vti diff --git a/Patterns/NumericalMethods/advection_LaxWendroff.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxWendroff.vti similarity index 100% rename from Patterns/NumericalMethods/advection_LaxWendroff.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxWendroff.vti diff --git a/Patterns/NumericalMethods/advection_LaxWendroff_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxWendroff_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_LaxWendroff_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_LaxWendroff_double.vti diff --git a/Patterns/NumericalMethods/advection_RungeKutta4.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_RungeKutta4.vti similarity index 100% rename from Patterns/NumericalMethods/advection_RungeKutta4.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_RungeKutta4.vti diff --git a/Patterns/NumericalMethods/advection_RungeKutta4_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_RungeKutta4_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_RungeKutta4_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_RungeKutta4_double.vti diff --git a/Patterns/NumericalMethods/advection_forwardEuler.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_forwardEuler.vti similarity index 100% rename from Patterns/NumericalMethods/advection_forwardEuler.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_forwardEuler.vti diff --git a/Patterns/NumericalMethods/advection_forwardEuler_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_forwardEuler_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_forwardEuler_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_forwardEuler_double.vti diff --git a/Patterns/NumericalMethods/advection_midpointMethod.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_midpointMethod.vti similarity index 100% rename from Patterns/NumericalMethods/advection_midpointMethod.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_midpointMethod.vti diff --git a/Patterns/NumericalMethods/advection_midpointMethod_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_midpointMethod_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_midpointMethod_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_midpointMethod_double.vti diff --git a/Patterns/NumericalMethods/advection_modifiedEuler.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_modifiedEuler.vti similarity index 100% rename from Patterns/NumericalMethods/advection_modifiedEuler.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_modifiedEuler.vti diff --git a/Patterns/NumericalMethods/advection_modifiedEuler_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_modifiedEuler_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_modifiedEuler_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_modifiedEuler_double.vti diff --git a/Patterns/NumericalMethods/advection_staggeredLeapfrog.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_staggeredLeapfrog.vti similarity index 100% rename from Patterns/NumericalMethods/advection_staggeredLeapfrog.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_staggeredLeapfrog.vti diff --git a/Patterns/NumericalMethods/advection_staggeredLeapfrog_double.vti b/Patterns/Experiments/TimHutton/NumericalMethods/advection_staggeredLeapfrog_double.vti similarity index 100% rename from Patterns/NumericalMethods/advection_staggeredLeapfrog_double.vti rename to Patterns/Experiments/TimHutton/NumericalMethods/advection_staggeredLeapfrog_double.vti From a9f3a2503e752103b8aa74b3a9583bc129b4490e Mon Sep 17 00:00:00 2001 From: Tim Hutton Date: Wed, 30 Dec 2020 11:05:18 +0000 Subject: [PATCH 20/20] Added note about how to view the two-slit experiment in 3D. Tweaked render settings. --- Patterns/Schrodinger1926/two_slit.vti | 31 ++++++++++++++------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/Patterns/Schrodinger1926/two_slit.vti b/Patterns/Schrodinger1926/two_slit.vti index 8fd28eb52..bdf228d5b 100644 --- a/Patterns/Schrodinger1926/two_slit.vti +++ b/Patterns/Schrodinger1926/two_slit.vti @@ -3,19 +3,18 @@ - The <a href="http://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation">Schrödinger equation</a> describes the change in the + The <a href="http://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation">Schrödinger equation</a> describes the change in the wave function for some quantum state. Here a wave packet travels through empty space until it reaches a high energy barrier with two slits. Some of the energy is reflected back but the rest makes a distinctive interference pattern. - + In our implementation, 'a' represents the real part of the complex number and 'b' the imaginary part. A fixed energy background is given by 'c'. The probability density is shown in 'd', with the energy barrier overlaid. - - By changing the values of parameter 'potential' and then using 'Generate Initial Pattern' from the 'Action' menu you can change the height - of the barrier. Or you can paint directly into 'c'. + + Change the dimensions to 64 x 64 x 64 to see the two-slit experiment in 3D. - + 0.001 4 @@ -24,11 +23,11 @@ delta_b = laplacian_a + c*a; d = a*a + b*b + c; - + - + - + @@ -46,7 +45,7 @@ - + @@ -56,7 +55,7 @@ - + @@ -73,7 +72,7 @@ - + @@ -81,9 +80,9 @@ - + - + @@ -95,8 +94,10 @@ + + - +