MAST
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
level_set_boundary_velocity.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 // MAST includes
24 #include "base/nonlinear_system.h"
25 
26 
28 MAST::FieldFunction<RealVectorX>("phi_vel"),
29 _dim (dim),
30 _phi (nullptr),
31 _mesh (nullptr),
32 _level_set_func (nullptr) {
33 
34 }
35 
37 
38 }
39 
40 
41 void
43  const MAST::MeshFieldFunction& phi) {
44 
45  libmesh_assert(!_phi);
46 
47  _phi = &phi;
48  _mesh = &sys.system().get_mesh();
49 }
50 
51 
52 
53 void
55 
56  _phi = nullptr;
57  _mesh = nullptr;
58  _level_set_func = nullptr;
59 }
60 
61 
62 void
64  const libMesh::Point& p,
65  const Real t,
66  RealVectorX& v) const {
67  this->velocity(f, p, t, v);
68 }
69 
70 
71 void
73  const libMesh::Point& p,
74  const Real t,
75  RealVectorX& v) const {
76 
77  libmesh_assert(_phi);
78 
80  val = RealVectorX::Zero(1),
81  grad = RealVectorX::Zero(_dim),
82  dval = RealVectorX::Zero(1);
83 
85  gradmat = RealMatrixX::Zero(1, _dim);
86 
87  // the velocity is identified using the level set function gradient
88  // and its sensitivity
89  (*_phi)(p, t, val);
90 
91  // since the function provides the velocity at the boundary, the
92  // level set value should be close to zero. This is commented out since
93  // for a coarse level set mesh used to define geometry on a fine analysis
94  // mesh, the boundary can exist on locations that are not necessarily phi=0.
95  //libmesh_assert_less_equal(std::fabs(val(0)), tol);
96  //if (std::fabs(val(0)) > tol)
97  // libMesh::out << "**** !!! level-set not zero at boundary: " << val(0) << std::endl;
98 
99  _phi->gradient(p, t, gradmat);
100  _phi->derivative(f, p, t, dval);
101 
102  // copy the matrix to the vector for further calculations
103  grad = gradmat.row(0);
104 
105  // at boundary, phi(x) = 0
106  // so, dphi/dp + grad(phi) . V = 0
107  // dphi/dp + grad(phi) . (-grad(phi)/|grad(phi)| Vn) = 0 [since V = -grad(phi)/|grad(phi)| Vn]
108  // dphi/dp -(grad(phi) . grad(phi))/|grad(phi)| Vn = 0
109  // Vn = (dphi/dp) |grad(phi)| / |grad(phi)|^2
110  // = (dphi/dp) / |grad(phi)|
111  // V = -grad(phi)/ |grad(phi)| Vn
112  // = -grad(phi)/ |grad(phi)| * (dphi/dp) / |grad(phi)|
113  // = -grad(phi) * (dphi/dp)/|grad(phi)|^2
114  //
115  v = -dval(0)/grad.squaredNorm()*grad;
116 }
117 
118 
119 void
122 
123  libmesh_assert(!_level_set_func);
124 
125  _level_set_func = &phi;
126 }
127 
128 
129 void
131 
132  _level_set_func = nullptr;
133 }
134 
135 
136 #include <fstream>
137 void
139 search_nearest_interface_point(const libMesh::Elem& e,
140  const unsigned int side,
141  const libMesh::Point& p,
142  const Real t,
143  RealVectorX& pt) const {
144 
145  libmesh_assert(_level_set_func);
146 
147  MAST::LevelSetIntersection intersection;
148 
149  const libMesh::Elem*
150  elem = e.neighbor_ptr(side);
151 
152  libmesh_assert(elem);
153 
154  intersection.init(*_level_set_func,
155  *elem, 0.,
156  _mesh->max_elem_id(),
157  _mesh->max_node_id());
158 
159  libMesh::Point p2;
160  intersection.get_nearest_intersection_point(p, p2);
161  pt(0) = p2(0); pt(1) = p2(1); pt(2) = p2(2);
162 // std::fstream o;
163 // o.open("pts.csv", std::fstream::app);
164 // o << p(0) << " , " << p(1) << " , 0 \n";
165 // o << pt(0) << " , " << pt(1) << " , 0 \n";
166 }
167 
168 
169 void
172  const libMesh::Elem& e,
173  const unsigned int side,
174  const libMesh::Point& p,
175  const Real t,
176  RealVectorX& v) const {
177 
178  libmesh_assert(_phi);
179 
180  // velocity at this interface point is given by the normal velocity
181  // So, we first find the point and then compute its velocity
182 
183  this->search_nearest_interface_point(e, side, p, t, v);
184 
185  libMesh::Point
186  pt(v(0), v(1), v(2));
187 
188  // now compute the velocity at this point.
189  v.setZero();
190  this->velocity(f, pt, t, v);
191 }
192 
193 
194 
195 
196 void
198 search_nearest_interface_point_old(const libMesh::Point& p,
199  const Real t,
200  const Real length,
201  RealVectorX& v,
202  bool allow_sub_search) const {
203 
204  libmesh_assert(_phi);
205 
207  phi = RealVectorX::Zero(1),
208  grad_phi = RealVectorX::Zero(_dim),
209  grad0 = RealVectorX::Zero(3),
210  grad1 = RealVectorX::Zero(3),
211  z = RealVectorX::Zero(3),
212  dval = RealVectorX::Zero(1),
213  p_ref = RealVectorX::Zero(3),
214  dv0 = RealVectorX::Zero(4),
215  dv1 = RealVectorX::Zero(4),
216  dx = RealVectorX::Zero(4),
217  gradL = RealVectorX::Zero(4);
218 
220  hess = RealMatrixX::Zero(3, 3),
221  gradmat = RealMatrixX::Zero(1, _dim),
222  coeffs = RealMatrixX::Zero(3, 3);
223 
224  libMesh::Point
225  p_opt = p;
226 
227  p_ref(0) = p(0);
228  p_ref(1) = p(1);
229  p_ref(2) = p(2);
230 
231  // The point is identified from a constrained optimization problem.
232  // p* = argmin { .5 * || p* - p ||^2: phi(p*) = 0 }
233  // This is formulated as a constrained minimization problem:
234  // (p*, l) = argmin { .5 * || p* - p ||^2 + l * phi(p*) },
235  // where the Lagrangian is defined as
236  // L (p*, l) = .5 * || p* - p ||^2 + l * phi(p*)
237  // = .5 * ((px* - px)^2 + (py* - py)^2 + (pz* - pz)^2) + l * phi(p*)
238  // The first order optimality for this provides
239  // grad ( L (p*, l)) = { (p*-p) + l * grad(phi(p*)); phi(p*) } = 0
240  // We solve this using Newton-Raphson with an approximate Hessian
241  //
242 
243  bool
244  if_cont = true;
245 
246  unsigned int
247  n_iters = 0,
248  max_iters = 80;
249 
250  Real
251  tol = 1.e-8,
252  L0 = 1.e12,
253  damp = 0.5,
254  L1 = 0.;
255 
256  // initialize the design point
257  dv0.topRows(3) = p_ref;
258  dv0(3) = 1.; // arbitrary value of Lagrange multiplier
259  dv1 = dv0;
260 
261  while (if_cont) {
262 
263  // compute the gradient of Lagrangian
264  (*_phi) (p_opt, t, phi);
265  _phi->gradient (p_opt, t, gradmat);
266 
267  // this defines the gradient of Lagrangian
268  gradL.topRows(3) = (dv0.topRows(3)-p_ref) + dv0(3) * gradmat.row(0).transpose();
269  gradL(3) = phi(0);
270 
271  // approximate hessian
272  coeffs.setZero(4,4);
273  for (unsigned int i=0; i<3; i++) {
274  coeffs(i,i) = 1.;
275  coeffs(3,i) = coeffs(i,3) = gradmat(0,i);
276  }
277  coeffs.topLeftCorner(3, 3) += hess;
278 
279  Eigen::FullPivLU<RealMatrixX> solver(coeffs);
280 
281  bool
282  continue_search = true;
283  dx = - solver.solve(gradL);
284  Real
285  factor = 1.;
286  unsigned int
287  n_it = 0,
288  max_it = 10;
289 
290 
291  // some searches may end up outside the mesh, and libMesh will through
292  // an exception. So, we catch the exception and take a smaller step.
293  while (continue_search) {
294  try {
295 
296  // update the design points and check for convergence
297  dv1 = dv0 + damp*factor * dx;
298  L1 = _evaluate_point_search_obj(p, t, dv1);
299  continue_search = false;
300  }
301  catch (...) {
302 
303  factor *= 0.5;
304  n_it ++;
305  if (n_it == max_it) {
306  // could not find a point inside the mesh. returning
307  // the reference point.
308  v.topRows(_dim) = dv0.topRows(_dim);
309  return;
310  }
311  else
312  continue_search = true;
313  }
314  }
315 
316  p_opt(0) = dv1(0); p_opt(1) = dv1(1); p_opt(2) = dv1(2);
317  L1 = _evaluate_point_search_obj(p, t, dv1);
318 
319  // update Jacobian
320  grad1 = gradmat.row(0);
321  z = (dv1.topRows(3)-dv0.topRows(3)) - hess * (grad1-grad0);
322  //hess += (z*z.transpose())/z.dot(grad1-grad0);
323  grad0 = grad1;
324  dv0 = dv1;
325  // reduce the step if the gradient is upated
326  /*if (z.norm())
327  damp = 0.3;
328  else
329  damp = 0.5;
330  */
331  if (n_iters == max_iters) {
332 
333  // instead, find the point closest to the latest point returned
334  // by the failed search. Do not allow another sub-search here
335  if (allow_sub_search) {
336  this->search_nearest_interface_point_old(p_opt, t, length, v, false);
337 
338  p_opt(0) = v(0); p_opt(1) = v(1); p_opt(2) = v(2);
339  dv0.topRows(3) = v;
340  (*_phi) (p_opt, t, phi);
341  }
342 
343  if_cont = false;
344  libMesh::Point dp = p_opt - p;
345  libMesh::out
346  << "Warning: nearest interface point search did not converge. Point found from sub-search. "
347  << std::endl;
348  libMesh::out
349  << " given pt: ("
350  << p(0) << " , " << p(1) << " , " << p(2)
351  << ") -> mapped pt: ("
352  << p_opt(0) << " , " << p_opt(1) << " , " << p_opt(2)
353  << "). phi = " << phi(0)
354  << " d/h = " << dp.norm()/length << std::endl;
355  }
356  if (std::fabs(gradL.norm()) <= tol) if_cont = false;
357  if (std::fabs((L1-L0)/L0) <= tol) if_cont = false;
358 
359  L0 = L1;
360  n_iters++;
361  }
362 
363  v.setZero();
364  v.topRows(_dim) = dv0.topRows(_dim);
365 }
366 
367 
368 
369 void
372  const libMesh::Point& p,
373  const Real t,
374  const Real length,
375  RealVectorX& v) const {
376 
377  libmesh_assert(_phi);
378 
379  // velocity at this interface point is given by the normal velocity
380  // So, we first find the point and then compute its velocity
381 
382  this->search_nearest_interface_point_old(p, t, length, v);
383 
384  libMesh::Point
385  pt(v(0), v(1), v(2));
386 
387  // now compute the velocity at this point.
388  v.setZero();
389  this->velocity(f, pt, t, v);
390 }
391 
392 
393 
394 
395 void
397  const Real t,
398  RealVectorX& n) const {
399 
400  libmesh_assert(_phi);
401 
403  gradmat = RealMatrixX::Zero(1, _dim);
404 
405  n.setZero();
406 
407  // the velocity is identified using the level set function gradient
408  // and its sensitivity
409  _phi->gradient(p, t, gradmat);
410  n.topRows(3) = -gradmat.row(0);
411 
412  n /= n.norm();
413 }
414 
415 
416 
417 void
419  const libMesh::Point& p,
420  const Real t,
421  RealVectorX& n) const {
422 
423  libmesh_assert(_phi);
424 
426  gradmat = RealMatrixX::Zero(1, _dim),
427  dgrad = RealMatrixX::Zero(1, _dim);
428 
429  n.setZero();
430 
431  // the velocity is identified using the level set function gradient
432  // and its sensitivity
433  _phi->gradient(p, t, gradmat);
434  _phi->perturbation_gradient(p, t, dgrad);
435 
437  v = gradmat.row(0),
438  dv = dgrad.row(0);
439 
440  n.topRows(3) = (-dv/v.norm() + v.dot(dv)/std::pow(v.dot(v),1.5) * v);
441 }
442 
443 
444 
445 Real
447 _evaluate_point_search_obj(const libMesh::Point& p,
448  const Real t,
449  const RealVectorX& dv) const {
450 
451  libMesh::Point p1(dv(0), dv(1), dv(2));
452  RealVectorX phi;
453  (*_phi)(p1, t, phi);
454  p1 -= p;
455  return .5 * p1.norm_sq() + dv(3)*phi(0);
456 }
457 
458 
virtual void perturbation_gradient(const libMesh::Point &p, const Real t, RealMatrixX &v) const
calculates the value of the function at the specified point, p, and time, t, and returns it in v...
virtual void gradient(const libMesh::Point &p, const Real t, RealMatrixX &g) const
calculates the gradient of value of the function at the specified point, p, and time, t, and returns it in g.
MAST::NonlinearSystem & system()
void init(MAST::SystemInitialization &sys, const MAST::MeshFieldFunction &phi)
virtual void derivative(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, RealVectorX &v) const
calculates the value of the derivative of function with respect to the function f at the specified po...
const MAST::MeshFieldFunction * _phi
This provides a wrapper FieldFunction compatible class that interpolates the solution using libMesh&#39;s...
virtual void derivative(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, RealVectorX &v) const
calculates the value of the function at the specified point, p, and time, t, and returns it in v...
void attach_level_set_function(const MAST::FieldFunction< Real > &phi)
attaches the level set function phi with this object.
LevelSetBoundaryVelocity(const unsigned int dim)
void search_nearest_interface_point_derivative_old(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, const Real length, RealVectorX &v) const
serches for a point pt in the vicinity of p on the level set interface, where level set function is z...
libMesh::Real Real
void normal_at_point(const libMesh::Point &p, const Real t, RealVectorX &n) const
void search_nearest_interface_point(const libMesh::Elem &e, const unsigned int side, const libMesh::Point &p, const Real t, RealVectorX &pt) const
void init(const MAST::FieldFunction< Real > &phi, const libMesh::Elem &e, const Real t, unsigned int max_elem_id, unsigned int max_node_id)
void velocity(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, RealVectorX &v) const
Matrix< Real, Dynamic, Dynamic > RealMatrixX
void search_nearest_interface_point_old(const libMesh::Point &p, const Real t, const Real length, RealVectorX &pt, bool allow_sub_search=true) const
serches for a point pt in the vicinity of p on the level set interface, where level set function is z...
This creates the base class for functions that have a saptial and temporal dependence, and provide sensitivity operations with respect to the functions and parameters.
const MAST::FieldFunction< Real > * _level_set_func
Real _evaluate_point_search_obj(const libMesh::Point &p, const Real t, const RealVectorX &dv) const
Matrix< Real, Dynamic, 1 > RealVectorX
void search_nearest_interface_point_derivative(const MAST::FunctionBase &f, const libMesh::Elem &e, const unsigned int side, const libMesh::Point &p, const Real t, RealVectorX &v) const
void get_nearest_intersection_point(const libMesh::Point &p, libMesh::Point &pt)
void clear_level_set_function()
clears the attached level set function
void normal_derivative_at_point(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, RealVectorX &n) const