MFEM  v4.6.0
Finite element discretization library
tmop_tools.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_TMOP_TOOLS_HPP
13 #define MFEM_TMOP_TOOLS_HPP
14 
15 #include "bilinearform.hpp"
16 #include "pbilinearform.hpp"
17 #include "tmop.hpp"
18 #include "gslib.hpp"
19 
20 namespace mfem
21 {
22 
23 // Performs the full remap advection loop.
25 {
26 private:
27  RK4Solver ode_solver;
28  Vector nodes0;
29  Vector field0;
30  const double dt_scale;
31  const AssemblyLevel al;
33 
34  void ComputeAtNewPositionScalar(const Vector &new_nodes, Vector &new_field);
35 public:
37  double timestep_scale = 0.5)
39  ode_solver(), nodes0(), field0(), dt_scale(timestep_scale), al(al) { }
40 
41  virtual void SetInitialField(const Vector &init_nodes,
42  const Vector &init_field);
43 
44  virtual void ComputeAtNewPosition(const Vector &new_nodes,
45  Vector &new_field,
46  int new_nodes_ordering = Ordering::byNODES);
47 
48  /// Set the memory type used for large memory allocations. This memory type
49  /// is used when constructing the AdvectorCGOper but currently only for the
50  /// parallel variant.
51  void SetMemoryType(MemoryType mt) { opt_mt = mt; }
52 };
53 
54 #ifdef MFEM_USE_GSLIB
56 {
57 private:
58  Vector nodes0;
59  GridFunction field0_gf;
60  FindPointsGSLIB *finder;
61 public:
62  InterpolatorFP() : finder(NULL) { }
63 
64  virtual void SetInitialField(const Vector &init_nodes,
65  const Vector &init_field);
66 
67  virtual void ComputeAtNewPosition(const Vector &new_nodes,
68  Vector &new_field,
69  int new_nodes_ordering = Ordering::byNODES);
70 
72  {
73  return finder;
74  }
75 
77  {
78  finder->FreeData();
79  delete finder;
80  }
81 };
82 #endif
83 
84 /// Performs a single remap advection step in serial.
86 {
87 protected:
88  const Vector &x0;
92  mutable BilinearForm M, K;
94 
95 public:
96  /** Here @a fes is the FESpace of the function that will be moved. Note
97  that Mult() moves the nodes of the mesh corresponding to @a fes. */
99  FiniteElementSpace &fes,
101 
102  virtual void Mult(const Vector &ind, Vector &di_dt) const;
103 };
104 
105 #ifdef MFEM_USE_MPI
106 /// Performs a single remap advection step in parallel.
108 {
109 protected:
110  const Vector &x0;
114  mutable ParBilinearForm M, K;
116 
117 public:
118  /** Here @a pfes is the ParFESpace of the function that will be moved. Note
119  that Mult() moves the nodes of the mesh corresponding to @a pfes.
120  @a mt is used to set the memory type of the integrators. */
121  ParAdvectorCGOper(const Vector &x_start, GridFunction &vel,
122  ParFiniteElementSpace &pfes,
125 
126  virtual void Mult(const Vector &ind, Vector &di_dt) const;
127 };
128 #endif
129 
131 {
132 protected:
133  // 0 - Newton, 1 - LBFGS.
135  bool parallel;
136 
137  // Line search step is rejected if min(detJ) <= min_detJ_threshold.
138  double min_detJ_threshold = 0.0;
139 
140  // Surface fitting variables.
141  mutable double surf_fit_err_avg_prvs = 10000.0;
143  mutable bool update_surf_fit_coeff = false;
144  double surf_fit_max_threshold = -1.0;
146  double surf_fit_scale_factor = 0.0;
147  mutable int adapt_inc_count = 0;
148  mutable int max_adapt_inc_count = 10;
149 
150  // Minimum determinant over the whole mesh. Used for mesh untangling.
151  double *min_det_ptr = nullptr;
152  // Flag to compute minimum determinant and maximum metric in ProcessNewState,
153  // which is required for TMOP_WorstCaseUntangleOptimizer_Metric.
154  mutable bool compute_metric_quantile_flag = true;
155 
156  // Quadrature points that are checked for negative Jacobians etc.
158  // These fields are relevant for mixed meshes.
161 
163 
165  {
166  if (IntegRules)
167  {
168  return IntegRules->Get(el.GetGeomType(), integ_order);
169  }
170  return ir;
171  }
172 
173  double ComputeMinDet(const Vector &x_loc,
174  const FiniteElementSpace &fes) const;
175 
176  double MinDetJpr_2D(const FiniteElementSpace*, const Vector&) const;
177  double MinDetJpr_3D(const FiniteElementSpace*, const Vector&) const;
178 
179  /** @name Methods for adaptive surface fitting weight. */
180  ///@{
181  /// Get the average and maximum surface fitting error at the marked nodes.
182  /// If there is more than 1 TMOP integrator, we get the maximum of the
183  /// average and maximum error over all integrators.
184  virtual void GetSurfaceFittingError(const Vector &x_loc,
185  double &err_avg, double &err_max) const;
186 
187  /// Update surface fitting weight as surf_fit_weight *= factor.
188  void UpdateSurfaceFittingWeight(double factor) const;
189 
190  /// Get the surface fitting weight for all the TMOP integrators.
191  void GetSurfaceFittingWeight(Array<double> &weights) const;
192  ///@}
193 
194 public:
195 #ifdef MFEM_USE_MPI
196  TMOPNewtonSolver(MPI_Comm comm, const IntegrationRule &irule, int type = 0)
197  : LBFGSSolver(comm), solver_type(type), parallel(true),
198  ir(irule), IntegRules(NULL), integ_order(-1) { }
199 #endif
200  TMOPNewtonSolver(const IntegrationRule &irule, int type = 0)
201  : LBFGSSolver(), solver_type(type), parallel(false),
202  ir(irule), IntegRules(NULL), integ_order(-1) { }
203 
204  /// Prescribe a set of integration rules; relevant for mixed meshes.
205  /** If called, this function has priority over the IntegrationRule given to
206  the constructor of the class. */
207  void SetIntegrationRules(IntegrationRules &irules, int order)
208  {
209  IntegRules = &irules;
210  integ_order = order;
211  }
212 
213  void SetMinDetPtr(double *md_ptr) { min_det_ptr = md_ptr; }
214 
215  /// Set the memory type for temporary memory allocations.
217 
218  /// Compute scaling factor for the node movement direction using line-search.
219  /// We impose constraints on TMOP energy, gradient, minimum Jacobian of
220  /// the mesh, and (optionally) on the surface fitting error.
221  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
222 
223  /// Update (i) discrete functions at new nodal positions, and
224  /// (ii) surface fitting weight.
225  virtual void ProcessNewState(const Vector &x) const;
226 
227  /** @name Methods for adaptive surface fitting weight. (Experimental) */
228  /// Enable/Disable adaptive surface fitting weight.
229  /// The weight is modified after each TMOPNewtonSolver iteration as:
230  /// w_{k+1} = w_{k} * @a surf_fit_scale_factor if relative change in
231  /// max surface fitting error < @a surf_fit_rel_change_threshold.
232  /// The solver terminates if the maximum surface fitting error does
233  /// not sufficiently decrease for @a max_adapt_inc_count consecutive
234  /// solver iterations or if the max error falls below @a surf_fit_max_threshold.
236  {
237  surf_fit_scale_factor = 10.0;
239  }
241  {
242  surf_fit_scale_factor = factor;
243  }
245  {
246  surf_fit_rel_change_threshold = threshold;
247  }
249  {
250  max_adapt_inc_count = count;
251  }
253  {
254  surf_fit_max_threshold = max_error;
255  }
256  void SetMinimumDeterminantThreshold(double threshold)
257  {
258  min_detJ_threshold = threshold;
259  }
260 
261  virtual void Mult(const Vector &b, Vector &x) const
262  {
263  if (solver_type == 0)
264  {
265  NewtonSolver::Mult(b, x);
266  }
267  else if (solver_type == 1)
268  {
269  LBFGSSolver::Mult(b, x);
270  }
271  else { MFEM_ABORT("Invalid type"); }
272  }
273 
274  virtual void SetSolver(Solver &solver)
275  {
276  if (solver_type == 0)
277  {
278  NewtonSolver::SetSolver(solver);
279  }
280  else if (solver_type == 1)
281  {
282  LBFGSSolver::SetSolver(solver);
283  }
284  else { MFEM_ABORT("Invalid type"); }
285  }
286  virtual void SetPreconditioner(Solver &pr) { SetSolver(pr); }
287 };
288 
289 void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm,
290  const TargetConstructor &tc, Mesh &pmesh,
291  char *title, int position);
292 #ifdef MFEM_USE_MPI
293 void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm,
294  const TargetConstructor &tc, ParMesh &pmesh,
295  char *title, int position);
296 #endif
297 
298 }
299 
300 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:233
void SetMinimumDeterminantThreshold(double threshold)
Definition: tmop_tools.hpp:256
VectorGridFunctionCoefficient u_coeff
Definition: tmop_tools.hpp:91
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
Definition: tmop_tools.cpp:879
virtual void GetSurfaceFittingError(const Vector &x_loc, double &err_avg, double &err_max) const
Definition: tmop_tools.cpp:665
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field, int new_nodes_ordering=Ordering::byNODES)
Definition: tmop_tools.cpp:370
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:207
const IntegrationRule & ir
Definition: tmop_tools.hpp:157
virtual void Mult(const Vector &ind, Vector &di_dt) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: tmop_tools.cpp:214
void GetSurfaceFittingWeight(Array< double > &weights) const
Get the surface fitting weight for all the TMOP integrators.
Definition: tmop_tools.cpp:635
double ComputeMinDet(const Vector &x_loc, const FiniteElementSpace &fes) const
Definition: tmop_tools.cpp:826
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:2007
Container class for integration rules.
Definition: intrules.hpp:412
double MinDetJpr_2D(const FiniteElementSpace *, const Vector &) const
Definition: tmop_pa_jp2.cpp:75
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:689
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
Definition: tmop_tools.cpp:907
void SetMinDetPtr(double *md_ptr)
Definition: tmop_tools.hpp:213
double surf_fit_rel_change_threshold
Definition: tmop_tools.hpp:145
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field, int new_nodes_ordering=Ordering::byNODES)
Definition: tmop_tools.cpp:29
void SetTerminationWithMaxSurfaceFittingError(double max_error)
Definition: tmop_tools.hpp:252
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
void SetAdaptiveSurfaceFittingRelativeChangeThreshold(double threshold)
Definition: tmop_tools.hpp:244
virtual void FreeData()
Definition: gslib.cpp:272
double b
Definition: lissajous.cpp:42
const AssemblyLevel al
Definition: tmop_tools.hpp:93
virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
Definition: tmop_tools.cpp:379
Performs a single remap advection step in serial.
Definition: tmop_tools.hpp:85
virtual void Mult(const Vector &ind, Vector &di_dt) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: tmop_tools.cpp:283
void EnableAdaptiveSurfaceFitting()
Definition: tmop_tools.hpp:235
TMOPNewtonSolver(const IntegrationRule &irule, int type=0)
Definition: tmop_tools.hpp:200
IntegrationRules * IntegRules
Definition: tmop_tools.hpp:159
const FindPointsGSLIB * GetFindPointsGSLIB() const
Definition: tmop_tools.hpp:71
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1819
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:163
virtual void ProcessNewState(const Vector &x) const
Definition: tmop_tools.cpp:706
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: tmop_tools.hpp:286
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
AdvectorCG(AssemblyLevel al=AssemblyLevel::LEGACY, double timestep_scale=0.5)
Definition: tmop_tools.hpp:36
const AssemblyLevel al
Definition: tmop_tools.hpp:115
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
TMOPNewtonSolver(MPI_Comm comm, const IntegrationRule &irule, int type=0)
Definition: tmop_tools.hpp:196
Performs a single remap advection step in parallel.
Definition: tmop_tools.hpp:107
double MinDetJpr_3D(const FiniteElementSpace *, const Vector &) const
Definition: tmop_pa_jp3.cpp:77
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:780
const IntegrationRule & GetIntegrationRule(const FiniteElement &el) const
Definition: tmop_tools.hpp:164
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void SetMemoryType(MemoryType mt)
Definition: tmop_tools.hpp:51
virtual void SetInitialField(const Vector &init_nodes, const Vector &init_field)
Definition: tmop_tools.cpp:22
VectorGridFunctionCoefficient u_coeff
Definition: tmop_tools.hpp:113
SerialAdvectorCGOper(const Vector &x_start, GridFunction &vel, FiniteElementSpace &fes, AssemblyLevel al=AssemblyLevel::LEGACY)
Definition: tmop_tools.cpp:193
void SetTempMemoryType(MemoryType mt)
Set the memory type for temporary memory allocations.
Definition: tmop_tools.hpp:216
Class for parallel bilinear form.
ParAdvectorCGOper(const Vector &x_start, GridFunction &vel, ParFiniteElementSpace &pfes, AssemblyLevel al=AssemblyLevel::LEGACY, MemoryType mt=MemoryType::DEFAULT)
Definition: tmop_tools.cpp:253
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points...
Definition: gslib.hpp:52
Vector data type.
Definition: vector.hpp:58
Vector coefficient defined by a vector GridFunction.
void SetMaxNumberofIncrementsForAdaptiveFitting(int count)
Definition: tmop_tools.hpp:248
Base class for solvers.
Definition: operator.hpp:682
void SetAdaptiveSurfaceFittingScalingFactor(double factor)
Definition: tmop_tools.hpp:240
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: tmop_tools.hpp:261
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: tmop_tools.hpp:274
void UpdateSurfaceFittingWeight(double factor) const
Update surface fitting weight as surf_fit_weight *= factor.
Definition: tmop_tools.cpp:610
virtual void SetInitialField(const Vector &init_nodes, const Vector &init_field)
Definition: tmop_tools.cpp:333