Line data Source code
1 : //
2 : // This is a general-purpose routine for specifying time-interpolated
3 : // quantities in an input file.
4 : // Interpolators obey the following syntax:
5 : //
6 : // .. code-block:: bash
7 : //
8 : // (t1,t2,..,tn;v1,v2,...,vn)
9 : //
10 : // where :code:`tn` are times and :code:`vn` are values at those times.
11 : //
12 : // :bdg-danger-line:`Note: do not include whitespace in an interpolator, this will mess with the parser.`
13 : //
14 : // Interpolators can usually be used in place of regular numbers, but only
15 : // when supported by the calling parse function.
16 : //
17 : // For instance, the `BC::Operator::Elastic::TensionTest` allows the user to specify
18 : // fixed displacement; for instance, to specify a value of 0.1:
19 : //
20 : // .. code-block:: bash
21 : //
22 : // bc.tenstion_test.disp = 0.1
23 : //
24 : // However, it may be desirable to change the value over time.
25 : // In this case, one may specify an interpolator string instead:
26 : //
27 : // .. code-block:: bash
28 : //
29 : // bc.tenstion_test.disp = (1.0,1.5,2.0;0.0,0.1,0.0)
30 : //
31 : // This will cause the displacement value to be zero for :code:`t<1.0`;
32 : // ramp to :code:`0.1` as :code:`1.0<t<2.0`;
33 : // ramp back to :code:`0.0` as :code:`2.0<t<3.0`;
34 : // and then remain at the constant value of :code:`3.0` for :code:`t>3.0`.
35 : //
36 :
37 : #ifndef NUMERIC_INTERPOLATOR_LINEAR_H_
38 : #define NUMERIC_INTERPOLATOR_LINEAR_H_
39 :
40 : #include <AMReX.H>
41 : #include <AMReX_MultiFab.H>
42 : #include "Set/Set.H"
43 : #include "Util/Util.H"
44 : #include "Numeric/Interpolator/Interpolator.H"
45 : #include "IO/ParmParse.H"
46 :
47 : namespace Numeric
48 : {
49 : namespace Interpolator
50 : {
51 : template <class T>
52 : class Linear : public Interpolator<T>
53 : {
54 : public:
55 482 : Linear() {};
56 694 : ~Linear() {};
57 496 : Linear(const Set::Scalar &a_val)
58 496 : {
59 496 : interval_points.push_back(0.0);
60 496 : data_points.push_back(a_val);
61 496 : }
62 104 : Linear(const std::string a_str)
63 104 : {
64 104 : define(a_str);
65 104 : }
66 2 : Linear(const std::vector<T> _data_points, const std::vector<Set::Scalar> _interval_points)
67 2 : {
68 2 : define(_data_points,_interval_points);
69 2 : };
70 :
71 106 : void define(const std::string a_str)
72 : {
73 212 : std::string str = a_str;
74 :
75 106 : if (Util::String::Contains(str,"(") || Util::String::Contains(str,":") || Util::String::Contains(str,")") || Util::String::Contains(str,","))
76 : {
77 36 : if (!Util::String::Contains(str,"(")) Util::Abort(INFO,"Mismatched parentheses while trying to parse ",str);
78 36 : if (!Util::String::Contains(str,")")) Util::Abort(INFO,"Mismatched parentheses while trying to parse ",str);
79 36 : if (!Util::String::Contains(str,":")) Util::Abort(INFO,"Missing : in interpolator while trying to parse ",str);
80 36 : Util::String::ReplaceAll(str,"(","");
81 36 : Util::String::ReplaceAll(str,")","");
82 :
83 72 : std::vector<std::string> splitstr = Util::String::Split(str,':');
84 72 : std::string str_time = splitstr[0];
85 72 : std::string str_val = splitstr[1];
86 :
87 72 : std::vector<std::string> str_time_arr = Util::String::Split(str_time,',');
88 72 : std::vector<std::string> str_val_arr = Util::String::Split(str_val,',');
89 :
90 36 : if (str_time_arr.size() != str_val_arr.size()) Util::Abort(INFO,"Mismatched number of time values vs array values while trying to parse ", str);
91 :
92 108 : for (unsigned int i = 0; i < str_time_arr.size(); i++)
93 : {
94 72 : interval_points.push_back(std::stod(str_time_arr[i]));
95 72 : if (i > 0) if (interval_points[i] < interval_points[i-1]) Util::Abort(INFO,"Time series must monotonically increase - error while trying to parse ", str);
96 72 : data_points.push_back(std::stod(str_val_arr[i]));
97 : }
98 : }
99 : else
100 : {
101 70 : interval_points.push_back(0.0);
102 70 : data_points.push_back(std::stod(a_str));
103 : }
104 106 : }
105 :
106 2 : void define(const std::vector<T> _data_points, const std::vector<Set::Scalar> _interval_points)
107 : {
108 2 : data_points = _data_points;
109 2 : interval_points = _interval_points;
110 2 : if(data_points.size() != interval_points.size())
111 0 : Util::Abort(INFO,"Data points and interval points have different sizes");
112 2 : };
113 :
114 77639820 : T operator() (const Set::Scalar point) const
115 : {
116 : // If there is only one point stored, return it.
117 77639820 : if (data_points.size() == 1) return data_points[0];
118 :
119 : // Do this if point is below interpolation range
120 9348 : if(point < interval_points[0])
121 : {
122 200 : return data_points[0];
123 : }
124 :
125 9148 : Set::Scalar interval = interval_points[1]-interval_points[0];
126 9148 : int start = 0;
127 9148 : if (uniform)
128 : {
129 9148 : start = (int)( (point - interval_points[0]) / interval);
130 9148 : start -= 1;
131 9148 : start = std::max(start,0);
132 9148 : start = std::min(start,(int)(interval_points.size()-2));
133 : }
134 :
135 : // Do this if point is within interpolation range
136 9742 : for (unsigned int i = start; i<interval_points.size(); i++)
137 : {
138 9544 : if (interval_points[i] <= point && point <= interval_points[i+1])
139 : {
140 8950 : return data_points[i] +
141 8950 : (point - interval_points[i]) * (data_points[i+1] - data_points[i]) /
142 8950 : (interval_points[i+1] - interval_points[i]);
143 :
144 : }
145 : }
146 :
147 : // Otherwise point must be outside interpolation range, so
148 : // return the value for the highest point
149 198 : return data_points[interval_points.size()-1];
150 : }
151 :
152 2 : static void Parse(Linear<T> & value, IO::ParmParse & pp)
153 : {
154 2 : std::string in;
155 2 : pp_query("str",in); // Interpolator string used when Parsed from queryclass.
156 2 : value.define(in);
157 2 : }
158 :
159 :
160 : static Linear<T> Read(std::string filename, int derivative = 0);
161 :
162 : protected:
163 : std::vector<T> data_points;
164 : std::vector<Set::Scalar> interval_points;
165 : bool uniform = true;
166 :
167 : public:
168 : template<class U> friend Linear<U> operator * (const Set::Scalar alpha, const Linear<U> &b);
169 : };
170 :
171 : template<class T>
172 : AMREX_FORCE_INLINE
173 : Linear<T> operator * (const Set::Scalar alpha, const Linear<T> &b)
174 : {
175 : Linear<T> ret;
176 : for (unsigned int i = 0; i < b.data_points.size(); i++)
177 : ret.data_points.push_back(alpha * b.data_points[i]);
178 : ret.interval_points = b.interval_points;
179 : ret.uniform = b.uniform;
180 : return ret;
181 : }
182 :
183 : template<>
184 : ALAMO_SINGLE_DEFINITION
185 0 : Linear<Set::Scalar> Linear<Set::Scalar>::Read(std::string filename, int derivative)
186 : {
187 0 : std::ifstream input;
188 0 : input.open(filename);
189 0 : std::string line;
190 0 : std::vector<Set::Scalar> theta, w;
191 : // parse input
192 0 : while (std::getline(input, line))
193 0 : { std::vector<std::string> dat = Util::String::Split(line);
194 0 : theta.push_back(std::stof(dat[0]));
195 0 : w.push_back(std::stof(dat[1]));
196 : }
197 0 : if (derivative == 0)
198 : {
199 0 : return Linear<Set::Scalar>(w,theta);
200 : }
201 0 : else if (derivative == 1)
202 : {
203 0 : std::vector<Set::Scalar> thetasmall, dw;
204 0 : for (unsigned int i = 1; i < theta.size() - 1; i++)
205 : {
206 0 : thetasmall.push_back(theta[i]);
207 0 : dw.push_back((w[i + 1] - w[i - 1]) / (theta[i + 1] - theta[i - 1]));
208 : }
209 0 : return Linear<Set::Scalar>(dw,thetasmall);
210 : }
211 0 : else if (derivative == 2)
212 : {
213 0 : std::vector<Set::Scalar> thetasmall, ddw;
214 0 : for (unsigned int i = 1; i < theta.size() - 1; i++)
215 : {
216 0 : thetasmall.push_back(theta[i]);
217 0 : ddw.push_back((w[i + 1] - 2.0 * w[i] + w[i - 1]) / ((theta[i + 1] - theta[i]) * (theta[i] - theta[i - 1])));
218 : }
219 0 : return Linear<Set::Scalar>(ddw,thetasmall);
220 : }
221 0 : else Util::Abort(INFO,"derivative must be 0,1,2 but got ",derivative);
222 0 : exit(-1);
223 : }
224 :
225 :
226 : }
227 : }
228 :
229 : #endif
|
