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