libMesh
biharmonic.C
Go to the documentation of this file.
1 // Libmesh includes
2 #include "libmesh/mesh_generation.h"
3 #include "libmesh/numeric_vector.h"
4 #include "libmesh/replicated_mesh.h"
5 
6 // Example includes
7 #include "biharmonic.h"
8 #include "biharmonic_jr.h"
9 
10 using namespace libMesh;
11 
12 // Constructor
14  EquationSystems(mesh),
15  _mesh(mesh)
16 {
17  // Retrieve parameters and set defaults
18  _verbose = false;
19  _growth = false;
20  _degenerate = false;
21  _cahn_hillard = false;
22  _netforce = false;
23 
24  if (on_command_line("--verbose"))
25  _verbose = true;
26  if (on_command_line("--growth"))
27  _growth = true;
28  if (on_command_line("--degenerate"))
29  _degenerate = true;
30  if (on_command_line("--cahn_hillard"))
31  _cahn_hillard = true;
32  if (on_command_line("--netforce"))
33  _netforce = true;
34 
35  _kappa = command_line_value("kappa", 1.0);
36 
37  // "type of energy (double well, double obstacle, logarithmic+double well, logarithmic+double obstacle)"
38  std::string energy = command_line_value("energy", std::string("double_well"));
39 
40  if (energy == "double_well")
42  else if (energy == "double_obstacle")
44  else if (energy == "log_double_well")
46  else if (energy == "log_double_obstacle")
48  else
49  libmesh_error_msg("Unknown energy type: " << energy);
50 
51  _tol = command_line_value("tol", 1.0e-8);
52  _theta = command_line_value("theta", .001);
53  _theta_c = command_line_value("theta_c", 1.0);
54 
55  // "order of log truncation (0=none, 2=quadratic, 3=cubic)"
56  _log_truncation = command_line_value("log_truncation", 2);
57 
58  if (!_log_truncation)
59  libMesh::out << "WARNING: no truncation is being used for the logarithmic free energy term.\nWARNING: division by zero possible!\n";
60 
61 
62  // Dimension
63  _dim = command_line_value("dim", 1);
64 
65  libmesh_assert_msg((_dim <= 3) && (_dim > 0), "Invalid mesh dimension");
66 
67  // Build the mesh
68  // Yes, it's better to make a coarse mesh and then refine it. We'll get to it later.
69  _N = command_line_value("N", 8);
70  libmesh_assert_msg(_N > 0, "Invalid mesh size");
71 
72  switch (_dim)
73  {
74  case 1:
76  break;
77  case 2:
78  MeshTools::Generation::build_square(_mesh, _N, _N, 0.0, 1.0, 0.0, 1.0, QUAD4);
79  break;
80  case 3:
81  MeshTools::Generation::build_cube(_mesh, _N, _N, _N, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, HEX8);
82  break;
83  default:
84  libmesh_assert_msg((_dim <= 3) && (_dim > 0), "Invalid mesh dimension");
85  break;
86  }
87 
88  // Determine the initial timestep size
89  _dt0 = command_line_value("dt", 1.0/(10*_kappa*_N*_N*_N*_N));
90  libmesh_assert_msg(_dt0>=0, "Negative initial timestep");
91 
92  _t0 = command_line_value("min_time", 0.0);
93  _t1 = command_line_value("max_time", _t0 + 50.0*_dt0);
94  libmesh_assert_msg(_t1 >= _t0, "Final time less than initial time");
95  _T = _t1 - _t0;
96 
97  _cnWeight = command_line_value("crank_nicholson_weight", 1.0);
98  libmesh_assert_msg(_cnWeight <= 1 && _cnWeight >= 0, "Crank-Nicholson weight must be between 0 and 1");
99 
100  // Initial state
102  std::string initialState = command_line_value("initial_state", std::string("strip"));
103 
104  if (initialState == std::string("ball"))
106  else if (initialState == std::string("strip"))
108  else if (initialState == std::string("rod"))
109  _initialState = ROD;
110  else
111  libmesh_error_msg("Unknown initial state: neither ball nor rod nor strip");
112 
113  std::vector<Real> icenter;
114  command_line_vector("initial_center", icenter);
115 
116  // Check that the point defining the center was in the right spatial dimension
117  if (icenter.size() > _dim)
118  libmesh_assert_msg(icenter.size() > _dim, "Invalid dimension for the initial state center of mass");
119 
120  // Pad
121  icenter.resize(3);
122  for (std::size_t i = icenter.size(); i < _dim; ++i)
123  icenter[i] = 0.5;
124 
125  for (unsigned int i = _dim; i < 3; ++i)
126  icenter[i] = 0.0;
127 
128  _initialCenter = Point(icenter[0], icenter[1], icenter[2]);
129  _initialWidth = command_line_value("initial_width", 0.125);
130 
131  // Build the main equation encapsulated in the JR (Jacobian-Residual or J(R) "jet of R") object
132  _jr = &(add_system<Biharmonic::JR>(std::string("Biharmonic::JR")));
133 
134  // Output options
135 #ifdef LIBMESH_HAVE_EXODUS_API
136  if (on_command_line("output_base"))
137  _ofile_base = command_line_value("output_base", std::string("bih"));
138 
139  else
140  {
141  switch(_dim)
142  {
143  case 1:
144  _ofile_base = std::string("bih.1");
145  break;
146  case 2:
147  _ofile_base = std::string("bih.2");
148  break;
149  case 3:
150  _ofile_base = std::string("bih.3");
151  break;
152  default:
153  _ofile_base = std::string("bih");
154  break;
155  }
156  }
157  _ofile = _ofile_base + ".e";
158  _exio.reset(new ExodusII_IO(_mesh));
159  _o_dt = command_line_value("output_dt", 0.0);
160 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
161 } // constructor
162 
163 
164 
166 {
167  libMesh::out << "Biharmonic parameters:\n";
168 
169  // Print verbosity status
170  if (_verbose)
171  libMesh::out << "verbose mode is on\n";
172  else
173  libMesh::out << "verbose mode is off\n";
174 
175  // Print parameters
176  libMesh::out << "mesh dimension = " << _dim << "\n";
177  libMesh::out << "initial linear mesh size = " << _N << "\n";
178  libMesh::out << "kappa = " << _kappa << "\n";
179  libMesh::out << "growth = " << (int)_growth << "\n";
180  libMesh::out << "degenerate = " << (int)_degenerate << "\n";
181  libMesh::out << "Cahn-Hillard = " << (int)_cahn_hillard << "\n";
182  libMesh::out << "netforce = " << (int)_netforce << "\n";
183  libMesh::out << "energy = " << _energy << "\n";
184  libMesh::out << "tol = " << _tol << "\n";
185  libMesh::out << "theta = " << _theta << "\n";
186  libMesh::out << "theta_c = " << _theta_c << "\n";
187  libMesh::out << "log truncation = " << _log_truncation << "\n";
188  libMesh::out << "initial timestep size = " << _dt0 << "\n";
189 
190  if (_initialState == STRIP)
191  libMesh::out << "initial state: strip\n";
192 
193  if (_initialState == ROD)
194  libMesh::out << "initial state: rod\n";
195 
196  if (_initialState == BALL)
197  libMesh::out << "initial state: ball\n";
198 
199  libMesh::out << "initial state center = " << _initialCenter(0) << "\n";
200  libMesh::out << "initial state width = " << _initialWidth << "\n";
201  libMesh::out << "initial time (min_time) = " << _t0 << "\n";
202  libMesh::out << "integration time = " << _T << "\n";
203  libMesh::out << "final time (max_time) = " << _t1 << "\n";
204  libMesh::out << "Crank-Nicholson weight = " << _cnWeight << "\n";
205  libMesh::out << "Output timestep = " << _o_dt << "\n";
206  libMesh::out << "Output filename base: " << _ofile_base << "\n";
207 }
208 
209 
210 
211 
213 {
214  if (_verbose)
215  libMesh::out << ">>> Initializing Biharmonic\n";
216 
217  _dt = 0;
218  _o_count = 0;
220 
221  if (_verbose)
222  libMesh::out << "<<< Initializing Biharmonic\n";
223 }
224 
225 
226 
227 
228 
229 void Biharmonic::step(const Real & dt_in)
230 {
231  // We need to update the old solution vector.
232  // The old solution vector will be the current solution vector from the
233  // previous time step. We use vector assignment. Only TransientSystems
234  // (and systems derived from them) contain old solutions.
235  if (dt_in < 0)
236  _dt = _dt0;
237  else
238  _dt = dt_in;
239 
240  *(_jr->old_local_solution) = *(_jr->current_local_solution);
241 
242  // this will localize the current solution, resulting in a
243  // current_local_solution with correct ghost values
244  _jr->solve();
245 }
246 
247 
248 
249 #ifdef LIBMESH_HAVE_EXODUS_API
250 void Biharmonic::output(int timestep,
251  const Real & t,
252  Real & o_t,
253  bool force)
254 {
255  if (!force && t - o_t < _o_dt)
256  return;
257 
258  ++_o_count;
259 
260  if (_verbose)
261  libMesh::out << "Writing state "
262  << timestep
263  << " at time "
264  << t
265  << " to file "
266  << _ofile
267  << "; output a total of "
268  << _o_count
269  << " states so far\n";
270 
271  _exio->write_timestep(_ofile, *this, timestep, t);
272 
273  if (!force)
274  o_t = t;
275 }
276 
277 #else
278 
279 // Avoid compiler warnings in non-standard build configurations.
280 void Biharmonic::output(int, const Real &, Real &, bool) {}
281 
282 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
283 
284 
285 
287 {
288  Real t = _t0, o_t = 0.0;
289  int timestep = 1;
290 
291  // Force-write the initial timestep
292  output(timestep, t, o_t, true);
293 
294  while (t < _t1)
295  {
296  ++timestep;
297 
298  // A pretty update message
299  if (_verbose)
300  libMesh::out << "Solving for state " << timestep << ", time " << t << "\n";
301 
302  // Move biharmonic one timestep forward
303  step();
304 
305  // Keep track of time and timestep
306  t += _dt;
307 
308  // Output
309  output(timestep, t, o_t);
310  } // while(t < _t1)
311 
312  // Force-write the final timestep
313  output(timestep, t, o_t, true);
314 }
void run()
Definition: biharmonic.C:286
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
std::string _ofile
Definition: biharmonic.h:106
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
bool _verbose
Definition: biharmonic.h:99
Biharmonic(ReplicatedMesh &mesh)
Constructor retrieves command-line options, setting defaults, if necessary.
Definition: biharmonic.C:13
Real _cnWeight
Definition: biharmonic.h:104
MeshBase & mesh
void init()
Initialize all the systems.
Definition: biharmonic.C:212
Real _kappa
Definition: biharmonic.h:94
Real _theta_c
Definition: biharmonic.h:94
void output(int timestep, const Real &t, Real &o_t, bool force=false)
Definition: biharmonic.C:250
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
The libMesh namespace provides an interface to certain functionality in the library.
Point _initialCenter
Definition: biharmonic.h:101
Real _tol
Definition: biharmonic.h:95
bool _netforce
Definition: biharmonic.h:96
int _o_count
Definition: biharmonic.h:109
void viewParameters()
Definition: biharmonic.C:165
bool _cahn_hillard
Definition: biharmonic.h:96
T command_line_value(const std::string &, T)
Definition: libmesh.C:932
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
FreeEnergyEnum _energy
Definition: biharmonic.h:97
int _log_truncation
Definition: biharmonic.h:98
unsigned int _dim
Definition: biharmonic.h:93
void step(const Real &dt=-1.0)
Definition: biharmonic.C:229
Real _initialWidth
Definition: biharmonic.h:102
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void command_line_vector(const std::string &, std::vector< T > &)
Definition: libmesh.C:978
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
OStreamProxy out
bool _growth
Definition: biharmonic.h:96
void build_line(UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 1D meshes.
Real _o_dt
Definition: biharmonic.h:108
unsigned int _N
Definition: biharmonic.h:93
UniquePtr< ExodusII_IO > _exio
Definition: biharmonic.h:107
ReplicatedMesh & _mesh
Definition: biharmonic.h:112
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
std::string _ofile_base
Definition: biharmonic.h:106
InitialStateEnum _initialState
Definition: biharmonic.h:100
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
Real _theta
Definition: biharmonic.h:94
bool _degenerate
Definition: biharmonic.h:96