fem_system_ex4.C 9.22 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
// The libMesh Finite Element Library.
// Copyright (C) 2002-2016 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

// <h1>FEMSystem Example 4 - Mixed-dimension heat transfer equation
// with FEMSystem</h1>
// \author Roy Stogner
// \date 2015
//
// This example shows how elements of multiple dimensions can be
// linked and computed upon simultaneously within the
// DifferentiableSystem class framework

// C++ includes
#include <iomanip>

// Basic include files
#include "libmesh/elem.h"
#include "libmesh/equation_systems.h"
#include "libmesh/error_vector.h"
#include "libmesh/exact_solution.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/getpot.h"
#include "libmesh/gmv_io.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/parsed_function.h"
#include "libmesh/uniform_refinement_estimator.h"

// The systems and solvers we may use
#include "heatsystem.h"
#include "libmesh/diff_solver.h"
#include "libmesh/euler_solver.h"
#include "libmesh/steady_solver.h"

// Bring in everything from the libMesh namespace
using namespace libMesh;

// The main program.
int main(int argc, char** argv) {
	// Initialize libMesh.
	LibMeshInit init(argc, argv);

#ifndef LIBMESH_ENABLE_AMR
	libmesh_example_requires(false, "--enable-amr");
#else

	// This doesn't converge with Eigen BICGSTAB for some reason...
	libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS,
	                         "--enable-petsc");

	// This doesn't converge without at least double precision
	libmesh_example_requires(sizeof(Real) > 4, "--disable-singleprecision");

	// Parse the input file
	GetPot infile("../fem_system_ex4.in");

	// Read in parameters from the input file
	const Real global_tolerance = infile("global_tolerance", 0.);
	const unsigned int nelem_target = infile("n_elements", 400);
	const Real deltat = infile("deltat", 0.005);
	const unsigned int coarsegridsize = infile("coarsegridsize", 20);
	const unsigned int coarserefinements = infile("coarserefinements", 0);
	const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10);
	const unsigned int dim = infile("dimension", 2);

	// Skip higher-dimensional examples on a lower-dimensional libMesh build
	libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");

	// We have only defined 2 and 3 dimensional problems
	libmesh_assert(dim == 2 || dim == 3);

	// Create a mesh, with dimension to be overridden later, distributed
	// across the default MPI communicator.
	Mesh mesh(init.comm());

	// And an object to refine it
	MeshRefinement mesh_refinement(mesh);
	mesh_refinement.coarsen_by_parents() = true;
	mesh_refinement.absolute_global_tolerance() = global_tolerance;
	mesh_refinement.nelem_target() = nelem_target;
	mesh_refinement.refine_fraction() = 0.3;
	mesh_refinement.coarsen_fraction() = 0.3;
	mesh_refinement.coarsen_threshold() = 0.1;

	// Use the MeshTools::Generation mesh generator to create a uniform
	// grid on the square or cube.  We crop the domain at y=2/3 to allow
	// for a homogeneous Neumann BC in our benchmark there.
	boundary_id_type bcid = 3; // +y in 3D

	mesh.read("../meshes/bridge.xda");
	{
		// Add boundary elements corresponding to the +y boundary of our
		// volume mesh
		std::set<boundary_id_type> bcids;
		bcids.insert(bcid);
		mesh.get_boundary_info().add_elements(bcids, mesh);
		mesh.prepare_for_use();
	}

	// To work around ExodusII file format limitations, we need elements
	// of different dimensionality to belong to different subdomains.
	// Our interior elements defaulted to subdomain id 0, so we'll set
	// boundary elements to subdomain 2.
	{
		const MeshBase::element_iterator end_el = mesh.elements_end();
		for (MeshBase::element_iterator el = mesh.elements_begin();
		     el != end_el; ++el) {
			Elem* elem = *el;
			if (elem->dim() < dim) elem->subdomain_id() = 2;
		}
	}

	mesh_refinement.uniformly_refine(coarserefinements);

	// Print information about the mesh to the screen.
	mesh.print_info();

	// Create an equation systems object.
	EquationSystems equation_systems(mesh);

	// Declare the system "Heat" and its variables.
	HeatSystem& system = equation_systems.add_system<HeatSystem>("Heat");

	// Solve this as a steady system
	system.time_solver = UniquePtr<TimeSolver>(new SteadySolver(system));

	// Initialize the system
	equation_systems.init();

	// Set the time stepping options
	system.deltat = deltat;

	// And the nonlinear solver options
	DiffSolver& solver = *(system.time_solver->diff_solver().get());
	solver.quiet = infile("solver_quiet", true);
	solver.verbose = !solver.quiet;
	solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
	solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
	solver.relative_residual_tolerance =
	    infile("relative_residual_tolerance", 0.0);
	solver.absolute_residual_tolerance =
	    infile("absolute_residual_tolerance", 0.0);

	// And the linear solver options
	solver.max_linear_iterations = infile("max_linear_iterations", 50000);
	solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3);

	// Print information about the system to the screen.
	equation_systems.print_info();

	// Adaptively solve the steady solution
	unsigned int a_step = max_adaptivesteps;
	for (; a_step != max_adaptivesteps; ++a_step) {
		system.solve();

		system.postprocess();

		ErrorVector error;

		UniquePtr<ErrorEstimator> error_estimator;

		// To solve to a tolerance in this problem we
		// need a better estimator than Kelly
		if (global_tolerance != 0.) {
			// We can't adapt to both a tolerance and a mesh
			// size at once
			libmesh_assert_equal_to(nelem_target, 0);

			UniformRefinementEstimator* u = new UniformRefinementEstimator;

			// The lid-driven cavity problem isn't in H1, so
			// lets estimate L2 error
			u->error_norm = L2;

			error_estimator.reset(u);
		} else {
			// If we aren't adapting to a tolerance we need a
			// target mesh size
			libmesh_assert_greater(nelem_target, 0);

			// Kelly is a lousy estimator to use for a problem
			// not in H1 - if we were doing more than a few
			// timesteps we'd need to turn off or limit the
			// maximum level of our adaptivity eventually
			error_estimator.reset(new KellyErrorEstimator);
		}

		error_estimator->estimate_error(system, error);

		// Print out status at each adaptive step.
		Real global_error = error.l2_norm();
		libMesh::out << "Adaptive step " << a_step << ": " << std::endl;

		if (global_tolerance != 0.)
			libMesh::out << "Global_error = " << global_error << std::endl;

		if (global_tolerance != 0.)
			libMesh::out << "Worst element error = " << error.maximum()
			             << ", mean = " << error.mean() << std::endl;

		if (global_tolerance != 0.) {
			// If we've reached our desired tolerance, we
			// don't need any more adaptive steps
			if (global_error < global_tolerance) break;
			mesh_refinement.flag_elements_by_error_tolerance(error);
		} else {
			// If flag_elements_by_nelem_target returns true, this
			// should be our last adaptive step.
			if (mesh_refinement.flag_elements_by_nelem_target(error)) {
				mesh_refinement.refine_and_coarsen_elements();
				equation_systems.reinit();
				a_step = max_adaptivesteps;
				break;
			}
		}

		// Carry out the adaptive mesh refinement/coarsening
		mesh_refinement.refine_and_coarsen_elements();
		equation_systems.reinit();

		libMesh::out << "Refined mesh to " << mesh.n_active_elem()
		             << " active elements and "
		             << equation_systems.n_active_dofs() << " active dofs."
		             << std::endl;
	}
	// Do one last solve if necessary
	if (a_step == max_adaptivesteps) {
		system.solve();

		system.postprocess();
	}

#ifdef LIBMESH_HAVE_EXODUS_API
	ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API

#ifdef LIBMESH_HAVE_GMV
	GMVIO(mesh).write_equation_systems("out.gmv", equation_systems);
#endif // #ifdef LIBMESH_HAVE_GMV

#ifdef LIBMESH_HAVE_FPARSER
	// Check that we got close to the analytic solution
	ExactSolution exact_sol(equation_systems);
	const std::string exact_str =
	    (dim == 2) ? "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)";
	ParsedFunction<Number> exact_func(exact_str);
	exact_sol.attach_exact_value(0, &exact_func);
	exact_sol.compute_error("Heat", "T");

	Number err = exact_sol.l2_error("Heat", "T");

	// Print out the error value
	libMesh::out << "L2-Error is: " << err << std::endl;

	libmesh_assert_less(libmesh_real(err), 2e-3);

#endif // #ifdef LIBMESH_HAVE_FPARSER

#endif // #ifndef LIBMESH_ENABLE_AMR

	// All done.
	return 0;
}