libMesh
miscellaneous_ex7.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Miscellaneous Example 7 - The PetscDMNonlinearSolver and
21 // Variational Inequality (VI) problems</h1>
22 // \author Dmitry Karpeyev
23 // \date 2012
24 //
25 // In this example, LibMesh interfaces directly with PETSc's
26 // variational inequality solver through PetscDMNonlinearSolver
27 // (available in PETSc-3.3.0 or above).
28 
29 // Example include files
30 #include "biharmonic.h"
31 
32 // Bring in everything from the libMesh namespace
33 using namespace libMesh;
34 
35 // Print usage information if requested on command line
36 void print_help(int argc, char ** argv);
37 
38 int main(int argc, char ** argv)
39 {
40  // Initialize libMesh.
41  LibMeshInit init (argc, argv);
42  if (on_command_line("--help"))
43  print_help(argc, argv);
44  else
45  {
46 #if !defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
47  libmesh_example_requires(false, "--enable-second");
48 #elif !defined(LIBMESH_ENABLE_PERIODIC)
49  libmesh_example_requires(false, "--enable-periodic");
50 #endif
51 
52  // This is a PETSc-specific solver
53  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
54 
55  const int dim = command_line_value("dim", 1);
56 
57  // Skip higher-dimensional examples on a lower-dimensional libMesh build
58  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
59 
60  // This example only works with ReplicatedMesh.
61  ReplicatedMesh mesh(init.comm());
62  Biharmonic biharmonic(mesh);
63  biharmonic.viewParameters();
64  biharmonic.init();
65  biharmonic.run();
66  }
67  return 0;
68 }
69 
70 
71 
72 
73 
74 void print_help(int, char ** argv)
75 {
76  libMesh::out << "This example solves the Cahn-Hillard equation with chemical potential f:\n"
77  << " u_t = \\div(M(u)\\grad f(u))\n"
78  << "Here we have\n"
79  << " u, -1 <= u <= 1 -- relative concentration (difference of two concentrations in a binary mixture) \n"
80  << " M, M >= 0 -- mobility of the mixture\n"
81  << " f = \\delta E/\\delta u -- variational derivative of the free energy functional E\n"
82  << " E = \\int[\\kappa/2 |\\grac u|^ + g(u)]\n"
83  << "where the gradient term is the interfacial energy density with \\kappa quantifying the energy of the interface,\n"
84  << "and g(u) is the bulk energy density\n"
85  << " g(u) = \\theta L(u) + \\theta_c W(u),\n"
86  << "L(u) is the (optional, in this model) logarithmic term corresponding to the entropy of the mixture:\n"
87  << " L(u) = (\\theta/2)[(1+u)\\ln((1+u)/2) + (1-u)\\ln((1-u)/2)],\n"
88  << "where \\theta is related to the Boltzmann factor k_B T - a proxy for the absolute temperature T.\n"
89  << "L can be optionally approximated ('truncated') using a quadratic or a cubic polynomial on [-1,1]\n"
90  << "W(u) is the (optional, in this model) potential promoting demixing. It can take the form of \n"
91  << "a 'double well' potential\n"
92  << " W(u) = \\theta_c (u^4/4 - u^2/2),\n"
93  << " or \n"
94  << "a 'double obstacle' potential\n"
95  << " W(u) = (\\theta_c/2)(1-u^2),\n"
96  << "where \\theta_c is the critical 'temperature'.\n"
97  << "Finally, mobility M can be constant of 'degenerate', by which we mean that M is varying with u and \n"
98  << "vanishing (degenerating) whenever u reaches critical values +1 or -1:\n"
99  << " M(u) = 1.0\n"
100  << " or\n"
101  << " M(u) = (1.0 - u^2)\n"
102  << "Degenerate mobility should generally be used only in conjunction with logarithmic free energy terms.\n\n"
103  << "The equation is solved on a periodic domain (in 1D, 2D or 3D)\n"
104  << "using a Galerkin formulation with C^1 elements approximating the H^2_{per} function space.\n\n"
105  << "\n-----------\n"
106  << "COMPILING: "
107  << "\n-----------\n"
108  << "Compile as follows (assuming libmesh has been built): \n"
109  << "METHOD=<method> make \n"
110  << "where <method> is dbg or opt.\n"
111  << "\n-----------\n"
112  << "HELP: "
113  << "\n-----------\n"
114  << "Print this help message:\n"
115  << argv[0] << " --help\n"
116  << "\n-----------\n"
117  << "RUNNING: "
118  << "\n-----------\n"
119  << "Run in serial with PETSc-3.4 and above as follows:\n"
120  << "\n"
121  << argv[0] << "\n"
122  << " [--verbose] dim=<1|2|3> N=<number_of_linear_elements> \n"
123  << " kappa=<kappa_value> growth=<yes|no> degenerate=<yes|no> [--cahn-hillard] \n"
124  << " [--netforce] energy=<double_well|double_obstacle|log_double_well|log_double_obstacle> log_truncation_order=<2|3> \n"
125  << " theta=<theta_value> theta_c=<theta_c_value> \n"
126  << " initial_state=<ball|rod|strip> initial_center='x [y [z]]' initial_width=<width> \n"
127  << " min_time=<initial_time> max_time=<final_time> dt=<timestep_size> crank_nicholson_weight=<between_0_and_1> \n"
128  << " output_base=<base_filename> output_dt=<output_timestep_size> [-snes_type vinewtonrsls | vinewtonssls] \n"
129  << "(with PETSc-3.3 use -snes_type virs | viss).\n"
130  << argv[0] << " --verbose \n"
131  << "is a pretty good start.\n"
132  << "\nModeling a 1D system with 2D or 3D (for a strip the second and third components of the center are immaterial):\n"
133  << argv[0]<< " --verbose dim=1 N=1024 initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6\n"
134  << argv[0]<< " --verbose dim=2 N=64 initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
135  << argv[0]<< " --verbose dim=3 N=32 initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
136  << "\n"
137  << "Modeling a 2D system with 3D (for a rod the third component of the center is immaterial) \n"
138  << argv[0]<< " --verbose dim=2 N=64 initial_state=rod initial_center='0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
139  << argv[0]<< " --verbose dim=3 N=32 initial_state=rod initial_center='0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
140  << "\n"
141  << "A 3D system with an initial ball in the center\n"
142  << argv[0] << " --verbose dim=3 N=32 initial_state=ball initial_center='0.5 0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
143  << "\n"
144  << "Add -snes_type vinewtonrsls to run the variational inequality version that ensures the solution is between -1.0 and 1.0 at all times.\n"
145  << "If using PETSc-3.3, run with -snes_type virs.\n"
146  << "\n"
147  << std::endl;
148 }
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
void print_help(int argc, char **argv)
unsigned int dim
MeshBase & mesh
The Biharmonic class encapsulates most of the data structures necessary to calculate the biharmonic r...
Definition: biharmonic.h:46
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
int main(int argc, char **argv)
SolverPackage default_solver_package()
Definition: libmesh.C:995
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.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
OStreamProxy out