fem_system_ex4.C 9.34 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
// 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"

Thomas Steinreiter's avatar
Thomas Steinreiter committed
51
52
// user defined literal for libmesh's Real
constexpr auto operator"" _r(long double n) { return libMesh::Real(n); }
53
54
55
56

// The main program.
int main(int argc, char** argv) {
	// Initialize libMesh.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
57
	libMesh::LibMeshInit init{argc, argv};
58
59
60
61
62
63

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

	// This doesn't converge with Eigen BICGSTAB for some reason...
Thomas Steinreiter's avatar
Thomas Steinreiter committed
64
65
	libmesh_example_requires(libMesh::default_solver_package() !=
	                             libMesh::EIGEN_SOLVERS,
66
67
68
	                         "--enable-petsc");

	// This doesn't converge without at least double precision
Thomas Steinreiter's avatar
Thomas Steinreiter committed
69
70
	libmesh_example_requires(sizeof(libMesh::Real) > 4,
	                         "--disable-singleprecision");
71
72

	// Parse the input file
Thomas Steinreiter's avatar
Thomas Steinreiter committed
73
	GetPot infile{"../fem_system_ex4.in"};
74
75

	// Read in parameters from the input file
Thomas Steinreiter's avatar
Thomas Steinreiter committed
76
77
78
79
80
81
82
	const auto& global_tolerance{infile("global_tolerance", 0._r)};
	const auto& nelem_target{infile("n_elements", 400u)};
	const auto& deltat{infile("deltat", 0.005_r)};
	const auto& coarsegridsize{infile("coarsegridsize", 20u)};
	const auto& coarserefinements{infile("coarserefinements", 0u)};
	const auto& max_adaptivesteps{infile("max_adaptivesteps", 10u)};
	const auto& dim{infile("dimension", 2u)};
83
84
85
86
87
88
89
90
91

	// 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.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
92
	libMesh::Mesh mesh{init.comm()};
93
94

	// And an object to refine it
Thomas Steinreiter's avatar
Thomas Steinreiter committed
95
	libMesh::MeshRefinement mesh_refinement{mesh};
96
97
98
99
100
101
102
103
104
105
	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.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
106
	libMesh::boundary_id_type bcid{3}; // +y in 3D
107
108

	mesh.read("../meshes/bridge.xda");
Thomas Steinreiter's avatar
Thomas Steinreiter committed
109
110
111
112
113

	// Add boundary elements corresponding to the +y boundary of our
	// volume mesh
	mesh.get_boundary_info().add_elements({bcid}, mesh);
	mesh.prepare_for_use();
114
115
116
117
118

	// 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.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
119
120
121
	std::for_each(mesh.elements_begin(), mesh.elements_end(), [&](auto elem) {
		if (elem->dim() < dim) elem->subdomain_id() = 2;
	});
122
123
124
125
126
127
128

	mesh_refinement.uniformly_refine(coarserefinements);

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

	// Create an equation systems object.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
129
	libMesh::EquationSystems equation_systems{mesh};
130
131

	// Declare the system "Heat" and its variables.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
132
	auto& system = equation_systems.add_system<HeatSystem>("Heat");
133
134

	// Solve this as a steady system
Thomas Steinreiter's avatar
Thomas Steinreiter committed
135
136
	system.time_solver = libMesh::UniquePtr<libMesh::TimeSolver>(
	    new libMesh::SteadySolver(system));
137
138
139
140
141
142
143
144

	// Initialize the system
	equation_systems.init();

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

	// And the nonlinear solver options
Thomas Steinreiter's avatar
Thomas Steinreiter committed
145
	auto& solver = *(system.time_solver->diff_solver().get());
146
147
	solver.quiet = infile("solver_quiet", true);
	solver.verbose = !solver.quiet;
Thomas Steinreiter's avatar
Thomas Steinreiter committed
148
149
	solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15u);
	solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3_r);
150
	solver.relative_residual_tolerance =
Thomas Steinreiter's avatar
Thomas Steinreiter committed
151
	    infile("relative_residual_tolerance", 0.0_r);
152
	solver.absolute_residual_tolerance =
Thomas Steinreiter's avatar
Thomas Steinreiter committed
153
	    infile("absolute_residual_tolerance", 0.0_r);
154
155

	// And the linear solver options
Thomas Steinreiter's avatar
Thomas Steinreiter committed
156
157
158
	solver.max_linear_iterations = infile("max_linear_iterations", 50000u);
	solver.initial_linear_tolerance =
	    infile("initial_linear_tolerance", 1.e-3_r);
159
160
161
162
163

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

	// Adaptively solve the steady solution
Thomas Steinreiter's avatar
Thomas Steinreiter committed
164
	auto a_step = max_adaptivesteps;
165
166
167
168
169
	for (; a_step != max_adaptivesteps; ++a_step) {
		system.solve();

		system.postprocess();

Thomas Steinreiter's avatar
Thomas Steinreiter committed
170
		libMesh::UniquePtr<libMesh::ErrorEstimator> error_estimator;
171
172
173
174
175
176
177
		// 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);

Thomas Steinreiter's avatar
Thomas Steinreiter committed
178
179
			libMesh::UniformRefinementEstimator* u{
			    new libMesh::UniformRefinementEstimator};
180
181
182

			// The lid-driven cavity problem isn't in H1, so
			// lets estimate L2 error
Thomas Steinreiter's avatar
Thomas Steinreiter committed
183
			u->error_norm = libMesh::L2;
184
185
186
187
188
189
190
191
192
193
194

			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
Thomas Steinreiter's avatar
Thomas Steinreiter committed
195
			error_estimator.reset(new libMesh::KellyErrorEstimator);
196
197
		}

Thomas Steinreiter's avatar
Thomas Steinreiter committed
198
		libMesh::ErrorVector error;
199
200
201
		error_estimator->estimate_error(system, error);

		// Print out status at each adaptive step.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
202
		libMesh::out << "Adaptive step " << a_step << ": " << '\n';
203

Thomas Steinreiter's avatar
Thomas Steinreiter committed
204
		const auto& global_error = error.l2_norm();
205
		if (global_tolerance != 0.)
Thomas Steinreiter's avatar
Thomas Steinreiter committed
206
			libMesh::out << "Global_error = " << global_error << '\n';
207
208
209

		if (global_tolerance != 0.)
			libMesh::out << "Worst element error = " << error.maximum()
Thomas Steinreiter's avatar
Thomas Steinreiter committed
210
			             << ", mean = " << error.mean() << '\n';
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233

		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 "
Thomas Steinreiter's avatar
Thomas Steinreiter committed
234
		             << equation_systems.n_active_dofs() << " active dofs.\n";
235
236
237
238
239
240
241
242
243
	}
	// Do one last solve if necessary
	if (a_step == max_adaptivesteps) {
		system.solve();

		system.postprocess();
	}

#ifdef LIBMESH_HAVE_EXODUS_API
Thomas Steinreiter's avatar
Thomas Steinreiter committed
244
245
	libMesh::ExodusII_IO(mesh).write_equation_systems("out.e",
	                                                  equation_systems);
246
247
248
#endif // #ifdef LIBMESH_HAVE_EXODUS_API

#ifdef LIBMESH_HAVE_GMV
Thomas Steinreiter's avatar
Thomas Steinreiter committed
249
	libMesh::GMVIO(mesh).write_equation_systems("out.gmv", equation_systems);
250
251
252
253
#endif // #ifdef LIBMESH_HAVE_GMV

#ifdef LIBMESH_HAVE_FPARSER
	// Check that we got close to the analytic solution
Thomas Steinreiter's avatar
Thomas Steinreiter committed
254
	libMesh::ExactSolution exact_sol(equation_systems);
255
256
	const std::string exact_str =
	    (dim == 2) ? "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)";
Thomas Steinreiter's avatar
Thomas Steinreiter committed
257
	libMesh::ParsedFunction<libMesh::Number> exact_func(exact_str);
258
259
260
	exact_sol.attach_exact_value(0, &exact_func);
	exact_sol.compute_error("Heat", "T");

Thomas Steinreiter's avatar
Thomas Steinreiter committed
261
	libMesh::Number err{exact_sol.l2_error("Heat", "T")};
262
263

	// Print out the error value
Thomas Steinreiter's avatar
Thomas Steinreiter committed
264
	libMesh::out << "L2-Error is: " << err << '\n';
265
266
267
268
269
270
271
272
273
274

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

#endif // #ifdef LIBMESH_HAVE_FPARSER

#endif // #ifndef LIBMESH_ENABLE_AMR

	// All done.
	return 0;
}