Add monotone spline options to utilties. #1403

MFraters · 2017-03-07T16:25:22Z

This add an algorithm which can produce monotone splines. The code was based on the javascript implementation on https://en.wikipedia.org/wiki/Monotone_cubic_interpolation.

gassmoeller

In general a nice functionality. I have some stylistic and performance comments (it looks like this function is called a lot in the S20RTS and Savani initial conditions). Also please add a test for the new functionality (dont know if you already have a model that could be used, but otherwise you could create a simple initial condition plugin that only assigns the depth-dependent temperature of an arbitrarily chosen spline).
Additionally, it looks like the code style regarding spaces of this class is very different from the rest of ASPECT. Would you mind going over the code once and unifying this? In particular adding spaces around equal signs and at reasonable places in equations.

gassmoeller · 2017-03-14T18:50:31Z

source/utilities.cc

      {
        assert(x.size()==y.size());
        m_x=x;
        m_y=y;
-        int   n=x.size();
-        for (int i=0; i<n-1; i++)
+        unsigned int   n=x.size();


Thanks for improving this 👍. This can also be const, and please remove the useless spaces between int and n. Also is it possible to rename n to something useful? Is it n_knots or something else?

gassmoeller · 2017-03-14T18:53:06Z

source/utilities.cc

+                  {
+                    const double dx=x[i+1]-x[i];
+                    const double dy=y[i+1]-y[i];
+                    dxs.push_back(dx);


Please modify this so that dys,dxs, and ms to already have the correct size (n-1) upon creation. push_back can be slow, because the whole vector needs to copied and new memory needs to be allocated.

gassmoeller · 2017-03-14T18:54:25Z

source/utilities.cc

+                  }
+
+                // get m_a parameter
+                m_c.resize(1);


Same here, m_c can be once resized to m_c.resize(dxs.size()+1)

Upon taking another look, is dxs.size()+1 equivalent to n? That would be even simpler.

gassmoeller · 2017-03-14T18:57:56Z

source/utilities.cc

+
+                    const double invDx = 1/dxs[i];
+                    const double common0 = c1 + m_c[i+1] - m0 - m0;
+                    m_b.push_back((m0 - c1 - common0) * invDx);


Same here. You can resize m_b and m_a once before the loop to a size of n

gassmoeller · 2017-03-14T18:58:16Z

source/utilities.cc

+                // Get b and c coefficients
+                for (unsigned int i = 0; i < m_c.size()-1; i++)
+                  {
+                    const double c1=m_c[i];


please correct spaces

MFraters · 2017-03-27T17:13:42Z

@gassmoeller I hope this adresses all your comments, or do you want me to restyle the whole spline code?
I will try do design a good test later.

MFraters · 2017-03-28T10:01:27Z

I have now changed the spacing in the whole class anyway. Is this the kind of spacing you where looking for? I have also added a test to test this whole spline class.

tjhei · 2017-03-28T15:33:55Z

include/aspect/utilities.h

@@ -280,7 +280,8 @@ namespace aspect
           */
          void set_points(const std::vector<double> &x,
                          const std::vector<double> &y,
-                          bool cubic_spline = true);
+                          bool cubic_spline = true,
+                          bool monotome_spline = false);


typo: "monotone"

tjhei · 2017-03-28T15:34:51Z

source/utilities.cc

+                  {
+                    dxs[i] = x[i+1]-x[i];
+                    dys[i] = y[i+1]-y[i];
+                    ms[i] = double(dys[i])/double(dxs[i]);


these are doubles already, so remove the cast

gassmoeller · 2017-03-28T16:14:57Z

Thanks! This looks ready to merge except for the two minor things Timo mentioned. Lets see what the tester says.

gassmoeller · 2017-03-28T16:15:05Z

/run-tests

MFraters · 2017-04-04T13:55:02Z

Any other comments? :)

bangerth · 2017-04-04T16:50:19Z

include/aspect/utilities.h

@@ -280,7 +280,8 @@ namespace aspect
           */
          void set_points(const std::vector<double> &x,
                          const std::vector<double> &y,
-                          bool cubic_spline = true);
+                          bool cubic_spline = true,
+                          bool monotone_spline = false);


Can you document this additional argument? Also maybe in the class documentation?

It's not clear to me whether you can have both cubic and monotone splines. Is this possible?

Oh, also mark both of the bool arguments as const.

bangerth · 2017-04-04T16:52:07Z

source/utilities.cc

+                /**
+                 * This monotone spline algorithm is based on the javascript version
+                 * at https://en.wikipedia.org/wiki/Monotone_cubic_interpolation. The
+                 * parameters from this algorithm prevents overshooting in the


prevents -> prevent

bangerth · 2017-04-04T16:53:42Z

tests/splines.cc

@@ -0,0 +1,55 @@
+/**
+ * This tests the Utilities::tk::spine class for the linear case,


spine -> spline

MFraters force-pushed the add_monotone_spline branch from 1d5a802 to 6ad9c89 Compare March 8, 2017 10:59

gassmoeller reviewed Mar 14, 2017

View reviewed changes

MFraters added 4 commits March 28, 2017 13:53

Add monotone spline options to utilties.

ab5bcff

Address comments

9a39f04

Change the spacing in tk class

874393f

Add test for splines class.

13f9768

MFraters force-pushed the add_monotone_spline branch from d07c967 to 13f9768 Compare March 28, 2017 11:55

tjhei reviewed Mar 28, 2017

View reviewed changes

Address comments Timo

67bf16f

bangerth reviewed Apr 4, 2017

View reviewed changes

MFraters added 2 commits April 5, 2017 18:31

Adress comments Wolfgang spline utilities

29108c4

Adress comments Wolfgang spline test

a983c4f

gassmoeller merged commit 26e088d into geodynamics:master Apr 14, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add monotone spline options to utilties. #1403

Add monotone spline options to utilties. #1403

MFraters commented Mar 7, 2017

gassmoeller left a comment

gassmoeller Mar 14, 2017

gassmoeller Mar 14, 2017

gassmoeller Mar 14, 2017

gassmoeller Mar 14, 2017

gassmoeller Mar 14, 2017

gassmoeller Mar 14, 2017

MFraters commented Mar 27, 2017

MFraters commented Mar 28, 2017

tjhei Mar 28, 2017

tjhei Mar 28, 2017

gassmoeller commented Mar 28, 2017

gassmoeller commented Mar 28, 2017

MFraters commented Apr 4, 2017

bangerth Apr 4, 2017

bangerth Apr 4, 2017

bangerth Apr 4, 2017

bangerth Apr 4, 2017

		@@ -0,0 +1,55 @@
		/**
		* This tests the Utilities::tk::spine class for the linear case,

Add monotone spline options to utilties. #1403

Add monotone spline options to utilties. #1403

Conversation

MFraters commented Mar 7, 2017

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MFraters commented Mar 27, 2017

MFraters commented Mar 28, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gassmoeller commented Mar 28, 2017

gassmoeller commented Mar 28, 2017

MFraters commented Apr 4, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment