<i>
This program was written for fun by David Neckels (NCAR) while working
at Sandia (on the Wyoming Express bus to and from Corrales each day).
-The main purpose was to better understand Euler flow.
+The main purpose was to better understand Euler flow.
The code solves the basic Euler equations of gas dynamics, by using a
fully implicit Newton iteration (inspired by Sandia's Aria code). The
code may be configured by an input deck to run different simulations
with the solution $\mathbf{w}=(\rho v_1,\ldots,\rho v_d,\rho,
E)^{\top}$ consisting of $\rho$ the fluid density, ${\mathbf v}=(v_1,\ldots v_d)^T$ the
flow velocity (and thus $\rho\mathbf v$ being the linear momentum
-density), and
+density), and
$E$ the energy density of the gas. We interpret the equations above as
$\partial_t \mathbf{w}_i + \nabla \cdot \mathbf{F}_i(\mathbf{w}) = \mathbf
G_i(\mathbf w)$, $i=1,\ldots,dim+2$.
\partial_t (\rho v_i) + \sum_{s=1}^d \frac{\partial(\rho v_i v_s +
\delta_{is} p)}{\partial x_s} &=& g_i \rho, \qquad i=1,\dots,d, \\
\partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\
- \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=&
+ \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=&
\rho \mathbf g \cdot \mathbf v.
@f}
These equations describe, respectively, the conservation of momentum,
On the boundary, we have to say what the outer trace $\mathbf{w}^-$ is.
Depending on the boundary condition, we prescribe either of the following:
<ul>
-<li> Inflow boundary: $\mathbf{w}^-$ is prescribed to be the desired value.
+<li> Inflow boundary: $\mathbf{w}^-$ is prescribed to be the desired value.
<li> Supersonic outflow boundary: $\mathbf{w}^- = \mathbf{w}^+$
<li> Subsonic outflow boundary: $\mathbf{w}^- = \mathbf{w}^+$ except that the energy variable
-is modified to support a prescribed pressure $p_o$, i.e.
+is modified to support a prescribed pressure $p_o$, i.e.
$\mathbf{w}^- =(\rho^+, \rho v_1^+, \dots, \rho v_d^+, p_o/(\gamma -1) + 0.5 \rho |\mathbf{v}^+|^2)$
<li> Reflective boundary: we set $\mathbf{w}^-$ so that $(\mathbf{v}^+ + \mathbf{v}^-) \cdot \mathbf{n} = 0$ and
$\rho^- = \rho^+,E^-=E^+$.
We use a time stepping scheme to substitute the time derivative in the
above equations. At each time step, our full discretization is thus
-that the residual applied to any test
+that the residual applied to any test
function $\mathbf z$ equals zero:
@f{eqnarray*}
-R(\mathbf{W}_{n+1})(\mathbf z) &=&
+R(\mathbf{W}_{n+1})(\mathbf z) &=&
\int_{\Omega} \left(\frac{\mathbf{w}_{n+1} - \mathbf{w}_n}{\delta t},
\mathbf{z}\right)
- \int_{\Omega} \left(\mathbf{F}(\tilde{\mathbf{w}}),
-\nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z})
+\nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z})
\\
&& \qquad
+
\int_{\partial \Omega} \left(\mathbf{H}(\tilde{\mathbf{w}}^+),
-\mathbf{w}^-(\tilde{\mathbf{w}}^+), \mathbf{n}), \mathbf{z}\right)
+\mathbf{w}^-(\tilde{\mathbf{w}}^+), \mathbf{n}), \mathbf{z}\right)
-
\int_{\partial \Omega} \left(\mathbf{G}(\mathbf{w}),
-\mathbf{z}\right)
+\mathbf{z}\right)
\\
& = & 0
@f}
With these choices, equating the residual to zero results in a
nonlinear system of equations which we solve the nonlinear system by a
-Newton iteration, i.e. by iterating
+Newton iteration, i.e. by iterating
@f{eqnarray*}
R'(\mathbf{W}^k,\delta \mathbf{W})(\mathbf z) & = & -
R(\mathbf{W}^{k})(\mathbf z) \qquad \qquad \forall \mathbg z\in V_h \\
the variables being functions of the independent variables. These
variables are used in an algorithm, and as the variables are used,
their sensitivities with respect to the degrees of freedom are
-continuously updated.
+continuously updated.
One can imagine that for the full Jacobian,
this could be prohibitively expensive: the number of independent variables are
@endcode
The output are the derivatives $\frac{\partial c(a,b)}{\partial a},
-\frac{\partial c(a,b)}{\partial b}$ of $c(a,b)=2a+\cos(ab)$ at $a=1,b=2$.
+\frac{\partial c(a,b)}{\partial b}$ of $c(a,b)=2a+\cos(ab)$ at $a=1,b=2$.
It should be noted that Sacado provides more auto-differentation capabilities than the small subset
used in this program. However, understanding the example above is
<h3> Adaptivity </h3>
The example uses an ad-hoc refinement indicator that shows some usefulness in shock-type problems, and
-in the downhill flow example included. We refine according to the squared gradient of the density.
+in the downhill flow example included. We refine according to the squared gradient of the density.
Hanging nodes are handled by computing the numerical flux across cells that are of differing
refinement levels, rather than using the ConstraintMatrix class as in
all other tutorial programs so far. In this way, the example combines
Further, we enforce a maximum number of refinement levels to keep refinement under check. It is the
author's experience that for adaptivity for a time dependent problem, refinement can easily lead the simulation to
a screeching halt, because of time step restrictions if the mesh
-becomes to fine in any part of the domain, if care is not taken. The amount of refinement is
+becomes too fine in any part of the domain, if care is not taken. The amount of refinement is
limited in the example by letting the user specify the
maximum level of refinement that will be present anywhere in the mesh. In this way, refinement
tends not to slow the simulation to a halt. This, of course, is purely a heuristic strategy, and
<h3>Input deck, initial and boundary conditions</h3>
We use an input file deck to drive the simulation. In this way, we can alter the boundary conditions
-and other important properties of the simulation without having to recompile. For more information on
+and other important properties of the simulation without having to recompile. For more information on
the format, look at the <a href="#Results">results section</a>, where we
describe an example input file in more detail.
components of the solution vectors and the equations.
<li>A namespace that deals with everything that has to do with run-time
- parameters.
+ parameters.
<li>The <code>ConservationLaw</code> class that deals with time stepping,
outer nonlinear and inner linear solves, assembling the linear systems, and
- the top-level logic that drives all this.
+ the top-level logic that drives all this.
</ul>
The reason for this approach is that it separates the various concerns in a
equations: One would simply re-implement the members of the
<code>EulerEquations</code> class for some other hyperbolic equation, or
augment the existing equations by additional ones (for example by advecting
-additional variables, or my adding chemistry, etc). Such modifications,
+additional variables, or by adding chemistry, etc). Such modifications,
however, would not affect the time stepping, or the nonlinear solvers if
correctly done, and consequently nothing in the <code>ConservationLaw</code>
would have to be modified.
// introduction. We group all this into a
// structure that defines everything that has
// to do with the flux. All members of this
- // structures are static, i.e. the structure
+ // structure are static, i.e. the structure
// has no actual state specified by instance
// member variables. The better way to do
// this, rather than a structure with all
struct EulerEquations
{
// @sect4{Component description}
-
+
// First a few variables that
// describe the various components of our
// solution vector in a generic way. This
.push_back (DataComponentInterpretation::component_is_scalar);
data_component_interpretation
.push_back (DataComponentInterpretation::component_is_scalar);
-
+
return data_component_interpretation;
}
-
-
+
+
// @sect4{Transformations between variables}
-
+
// Next, we define the gas
// constant. We will set it to 1.4
// in its definition immediately
return ((gas_gamma-1.0) *
(*(W.begin() + energy_component) -
compute_kinetic_energy<number>(W)));
- }
+ }
+
+ // @sect4{EulerEquations::compute_flux_matrix}
- // @sect4{EulerEquations::compute_flux_matrix}
-
// We define the flux function
// $F(W)$ as one large matrix.
// Each row of this matrix
// matrix that correspond to the
// momentum terms:
const number pressure = compute_pressure<number> (W);
-
+
for (unsigned int d=0; d<dim; ++d)
{
for (unsigned int e=0; e<dim; ++e)
= W[first_momentum_component+d] *
W[first_momentum_component+e] /
W[density_component];
-
+
flux[first_momentum_component+d][d] += pressure;
}
-
+
// Then the terms for the
// density (i.e. mass
// conservation), and,
// lastly, conservation of
// energy:
for (unsigned int d=0; d<dim; ++d)
- flux[density_component][d] = W[first_momentum_component+d];
+ flux[density_component][d] = W[first_momentum_component+d];
for (unsigned int d=0; d<dim; ++d)
flux[energy_component][d] = W[first_momentum_component+d] /
{
Sacado::Fad::DFad<double> iflux[n_components][dim];
Sacado::Fad::DFad<double> oflux[n_components][dim];
-
+
compute_flux_matrix (Wplus, iflux);
compute_flux_matrix (Wminus, oflux);
-
+
for (unsigned int di=0; di<n_components; ++di)
{
normal_flux[di] = 0;
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int d=0; d<dim; ++d)
normal_flux[di] += 0.5*(iflux[di][d] + oflux[di][d]) * normal[d];
-
+
normal_flux[di] += 0.5*alpha*(Wplus[di] - Wminus[di]);
}
}
{
Wminus[c] = Wplus[c];
break;
- }
-
+ }
+
// Prescribed pressure boundary
// conditions are a bit more
// complicated by the fact that
Wminus[c] = boundary_values(c) / (gas_gamma-1.0) +
kinetic_energy;
-
+
break;
}
default:
Assert (false, ExcNotImplemented());
- }
+ }
}
compute_refinement_indicators (const DoFHandler<dim> &dof_handler,
const Mapping<dim> &mapping,
const Vector<double> &solution,
- Vector<double> &refinement_indicators)
- {
+ Vector<double> &refinement_indicators)
+ {
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
std::vector<unsigned int> dofs (dofs_per_cell);
std::vector<std::vector<Tensor<1,dim> > >
dU (1, std::vector<Tensor<1,dim> >(n_components));
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
= std::log(1+
std::sqrt(dU[0][density_component] *
dU[0][density_component]));
- }
+ }
}
-
-
-
+
+
+
// @sect4{EulerEquations::Postprocessor}
-
+
// Finally, we declare a class that
// implements a postprocessing of data
// components. The problem this class
// possibility to generate schlieren
// plots. Schlieren plots are a way to
// visualize shocks and other sharp
- // interfaces. The word "schlieren" a
+ // interfaces. The word "schlieren" is a
// German word that may be translated as
// "striae" -- it may be simpler to
// explain it by an example, however:
// schlieren is what you see when you,
// for example, pour highly concentrated
// alcohol, or a transparent saline
- // solution into water; the two have the
+ // solution, into water; the two have the
// same color, but they have different
// refractive indices and so before they
// are fully mixed light goes through the
// brightness variations if you look at
// it. That's "schlieren". A similar
// effect happens in compressible flow
- // due because the refractive index
+ // because the refractive index
// depends on the pressure (and therefore
// the density) of the gas.
//
{
public:
Postprocessor (const bool do_schlieren_plot);
-
+
virtual
void
compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
std::vector<Vector<double> > &computed_quantities) const;
virtual std::vector<std::string> get_names () const;
-
+
virtual
std::vector<DataComponentInterpretation::DataComponentInterpretation>
get_data_component_interpretation () const;
-
+
virtual UpdateFlags get_needed_update_flags () const;
virtual unsigned int n_output_variables() const;
const bool do_schlieren_plot;
};
};
-
+
template <int dim>
const double EulerEquations<dim>::gas_gamma = 1.4;
EulerEquations<dim>::Postprocessor::
get_needed_update_flags () const
{
- if (do_schlieren_plot == true)
+ if (do_schlieren_plot == true)
return update_values | update_gradients;
else
return update_values;
// perturbation that is added to the
// diagonal before forming the prec,
// and RTOL is a scaling factor $rtol
- // >= 1$.
+ // \geq 1$.
// - ilut_drop: The ILUT will
// drop any values that
// have magnitude less than this value.
// ParameterHandler::declare_entry
// call in
// <code>declare_parameters()</code>.
- struct Solver
+ struct Solver
{
enum SolverType { gmres, direct };
SolverType solver;
-
+
enum OutputType { quiet, verbose };
OutputType output;
}
prm.leave_subsection();
}
-
-
-
+
+
+
void Solver::parse_parameters (ParameterHandler &prm)
{
prm.enter_subsection("linear solver");
output = verbose;
if (op == "quiet")
output = quiet;
-
+
const std::string sv = prm.get("method");
- if (sv == "direct")
+ if (sv == "direct")
solver = direct;
else if (sv == "gmres")
solver = gmres;
ilut_rtol = prm.get_double("ilut relative tolerance");
ilut_drop = prm.get_double("ilut drop tolerance");
}
- prm.leave_subsection();
+ prm.leave_subsection();
}
-
-
+
+
// @sect4{Parameters::Refinement}
//
// Similarly, here are a few parameters
}
prm.leave_subsection();
}
-
+
void Refinement::parse_parameters (ParameterHandler &prm)
{
}
prm.leave_subsection();
}
-
+
// @sect4{Parameters::Flux}
// dependent value. In the latter case, it
// is chosen as $\frac{h}{2\delta T}$ with
// $h$ the diameter of the face to which
- // the flux is applied, and $\delta T$
+ // the flux is applied, and $\delta T$
// the current time step.
struct Flux
{
enum StabilizationKind { constant, mesh_dependent };
StabilizationKind stabilization_kind;
-
+
double stabilization_value;
static void declare_parameters (ParameterHandler &prm);
Patterns::Double(),
"alpha stabilization");
}
- prm.leave_subsection();
+ prm.leave_subsection();
}
-
-
+
+
void Flux::parse_parameters (ParameterHandler &prm)
{
prm.enter_subsection("flux");
{
const std::string stab = prm.get("stab");
- if (stab == "constant")
+ if (stab == "constant")
stabilization_kind = constant;
- else if (stab == "mesh")
+ else if (stab == "mesh")
stabilization_kind = mesh_dependent;
else
AssertThrow (false, ExcNotImplemented());
-
+
stabilization_value = prm.get_double("stab value");
}
prm.leave_subsection();
void Output::declare_parameters (ParameterHandler &prm)
{
prm.enter_subsection("output");
- {
+ {
prm.declare_entry("schlieren plot", "true",
Patterns::Bool (),
"Whether or not to produce schlieren plots");
}
prm.leave_subsection();
}
-
+
void Output::parse_parameters (ParameterHandler &prm)
{
typename EulerEquations<dim>::BoundaryKind
kind[EulerEquations<dim>::n_components];
-
+
FunctionParser<dim> values;
BoundaryConditions ();
};
-
-
+
+
AllParameters ();
-
+
double diffusion_power;
double time_step, final_time;
FunctionParser<dim> initial_conditions;
BoundaryConditions boundary_conditions[max_n_boundaries];
-
+
static void declare_parameters (ParameterHandler &prm);
void parse_parameters (ParameterHandler &prm);
};
:
initial_conditions (EulerEquations<dim>::n_components)
{}
-
+
template <int dim>
void
"outflow",
Patterns::Selection("inflow|outflow|pressure"),
"<inflow|outflow|pressure>");
-
+
prm.declare_entry("w_" + Utilities::int_to_string(di) +
" value", "0.0",
Patterns::Anything(),
}
else
is_stationary = false;
-
+
final_time = prm.get_double("final time");
theta = prm.get_double("theta scheme value");
}
{
std::vector<std::string>
expressions(EulerEquations<dim>::n_components, "0.0");
-
+
const bool no_penetration = prm.get_bool("no penetration");
for (unsigned int di=0; di<EulerEquations<dim>::n_components; ++di)
AssertThrow (false, ExcNotImplemented());
expressions[di] = prm.get("w_" + Utilities::int_to_string(di) +
- " value");
+ " value");
}
boundary_conditions[boundary_id].values
Parameters::Refinement::parse_parameters (prm);
Parameters::Flux::parse_parameters (prm);
Parameters::Output::parse_parameters (prm);
- }
+ }
}
-
-
+
+
// @sect3{Conservation law class}
public:
ConservationLaw (const char *input_filename);
void run ();
-
+
private:
void setup_system ();
// not of sufficiently high
// order.
Triangulation<dim> triangulation;
- const MappingQ1<dim> mapping;
-
+ const MappingQ1<dim> mapping;
+
const FESystem<dim> fe;
DoFHandler<dim> dof_handler;
const QGauss<dim> quadrature;
const QGauss<dim-1> face_quadrature;
-
+
// Next come a number of data
// vectors that correspond to the
// solution of the previous time
Epetra_SerialComm communicator;
std::auto_ptr<Epetra_Map> Map;
std::auto_ptr<Epetra_CrsMatrix> Matrix;
-
+
Parameters::AllParameters<dim> parameters;
ConditionalOStream verbose_cout;
};
dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
sparsity_pattern.compress();
-
+
std::vector<int> row_lengths (dof_handler.n_dofs());
for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
row_lengths[i] = sparsity_pattern.row_length (i);
const unsigned int max_nonzero_entries
= *std::max_element (row_lengths.begin(), row_lengths.end());
-
+
std::vector<double> values(max_nonzero_entries, 0);
std::vector<int> row_indices(max_nonzero_entries);
-
+
for (unsigned int row=0; row<dof_handler.n_dofs(); ++row)
{
row_indices.resize (row_lengths[row], 0);
values.resize (row_lengths[row], 0.);
-
+
for (int i=0; i<row_lengths[row]; ++i)
row_indices[i] = sparsity_pattern.column_number (row, i);
-
+
Matrix->InsertGlobalValues(row, row_lengths[row],
&values[0], &row_indices[0]);
}
| update_JxW_values
| update_normal_vectors,
neighbor_face_update_flags = update_values;
-
+
FEValues<dim> fe_v (mapping, fe, quadrature,
update_flags);
FEFaceValues<dim> fe_v_face (mapping, fe, face_quadrature,
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- for (; cell!=endc; ++cell)
+ for (; cell!=endc; ++cell)
{
fe_v.reinit (cell);
cell->get_dof_indices (dof_indices);
{
const unsigned int neighbor2=
cell->neighbor_of_neighbor(face_no);
-
+
for (unsigned int subface_no=0;
subface_no<GeometryInfo<dim>::subfaces_per_face;
++subface_no)
dof_indices_neighbor,
false,
numbers::invalid_unsigned_int,
- neighbor_child->diameter());
+ neighbor_child->diameter());
}
}
fe_v_subface_neighbor.reinit (neighbor,
neighbor_face_no,
neighbor_subface_no);
-
+
assemble_face_term (face_no, fe_v_face,
fe_v_subface_neighbor,
dof_indices,
cell->face(face_no)->diameter());
}
}
- }
+ }
// After all this assembling, notify the
// Trilinos matrix object that the matrix
void
ConservationLaw<dim>::
assemble_cell_term (const FEValues<dim> &fe_v,
- const std::vector<unsigned int> &dof_indices)
+ const std::vector<unsigned int> &dof_indices)
{
const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
const unsigned int n_q_points = fe_v.n_quadrature_points;
Table<2,Sacado::Fad::DFad<double> >
W_theta (n_q_points, EulerEquations<dim>::n_components);
-
+
Table<3,Sacado::Fad::DFad<double> >
grad_W (n_q_points, EulerEquations<dim>::n_components, dim);
std::vector<double> residual_derivatives (dofs_per_cell);
-
+
// Next, we have to define the independent
// variables that we will try to determine
// by solving a Newton step. These
std::vector<Sacado::Fad::DFad<double> > independent_local_dof_values(dofs_per_cell);
for (unsigned int i=0; i<dofs_per_cell; ++i)
independent_local_dof_values[i] = current_solution(dof_indices[i]);
-
+
// The next step incorporates all the
// magic: we declare a subset of the
// autodifferentiation variables as
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int c = fe_v.get_fe().system_to_component_index(i).first;
-
+
W[q][c] += independent_local_dof_values[i] *
fe_v.shape_value_component(i, q, c);
W_old[q][c] += old_solution(dof_indices[i]) *
(1-parameters.theta) *
old_solution(dof_indices[i])) *
fe_v.shape_value_component(i, q, c);
-
+
for (unsigned int d = 0; d < dim; d++)
grad_W[q][c][d] += independent_local_dof_values[i] *
fe_v.shape_grad_component(i, q, c)[d];
typedef Sacado::Fad::DFad<double> ForcingVector[EulerEquations<dim>::n_components];
ForcingVector *forcing = new ForcingVector[n_q_points];
-
+
for (unsigned int q=0; q<n_q_points; ++q)
{
EulerEquations<dim>::compute_flux_matrix (W_theta[q], flux[q]);
EulerEquations<dim>::compute_forcing_vector (W_theta[q], forcing[q]);
}
-
+
// We now have all of the pieces in place,
// so perform the assembly. We have an
// entries. Then, when we sum into the
// <code>right_hand_side</code> vector,
// we negate this residual.
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+ for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
{
Sacado::Fad::DFad<double> F_i = 0;
const unsigned int
component_i = fe_v.get_fe().system_to_component_index(i).first;
- // The residual for each row (i) will be accumulating
+ // The residual for each row (i) will be accumulating
// into this fad variable. At the end of the assembly
// for this row, we will query for the sensitivities
// to this variable and add them into the Jacobian.
fe_v.shape_value_component(i, point, component_i) *
fe_v.JxW(point);
- for (unsigned int d=0; d<dim; d++)
+ for (unsigned int d=0; d<dim; d++)
F_i -= flux[point][component_i][d] *
fe_v.shape_grad_component(i, point, component_i)[d] *
fe_v.JxW(point);
const std::vector<unsigned int> &dof_indices_neighbor,
const bool external_face,
const unsigned int boundary_id,
- const double face_diameter)
+ const double face_diameter)
{
const unsigned int n_q_points = fe_v.n_quadrature_points;
const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
const unsigned int n_independent_variables = (external_face == false ?
2 * dofs_per_cell :
dofs_per_cell);
-
+
for (unsigned int i = 0; i < dofs_per_cell; i++)
{
independent_local_dof_values[i] = current_solution(dof_indices[i]);
// of $W^-$ will be either functions of
// $W^+$, or they will be prescribed,
// depending on the kind of boundary
- // condition imposed here.
+ // condition imposed here.
//
// To start the evaluation, let us ensure
// that the boundary id specified for this
Wminus[q]);
}
-
+
// Now that we have $\mathbf w^+$ and
// $\mathbf w^-$, we can go about computing
// the numerical flux function $\mathbf
{
const unsigned int
component_i = fe_v.get_fe().system_to_component_index(i).first;
-
+
F_i += normal_fluxes[point][component_i] *
fe_v.shape_value_component(i, point, component_i) *
fe_v.JxW(point);
- }
+ }
for (unsigned int k=0; k<dofs_per_cell; ++k)
residual_derivatives[k] = F_i.fastAccessDx(k);
reinterpret_cast<int*>(
const_cast<unsigned int*>(
&dof_indices[0])));
-
+
if (external_face == false)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
const_cast<unsigned int*>(
&dof_indices_neighbor[0])));
}
-
+
right_hand_side(dof_indices[i]) -= F_i.val();
}
template <int dim>
std::pair<unsigned int, double>
-ConservationLaw<dim>::solve (Vector<double> &newton_update)
+ConservationLaw<dim>::solve (Vector<double> &newton_update)
{
Epetra_Vector x(View, *Map, newton_update.begin());
Epetra_Vector b(View, *Map, right_hand_side.begin());
-
+
switch (parameters.solver)
{
// If the parameter file specified that
solver.TrueResidual());
}
}
-
+
Assert (false, ExcNotImplemented());
return std::make_pair<unsigned int, double> (0,0);
}
if ((cell->level() < parameters.shock_levels) &&
(std::fabs(refinement_indicators(cell_no)) > parameters.shock_val))
cell->set_refine_flag();
- else
+ else
if ((cell->level() > 0) &&
(std::fabs(refinement_indicators(cell_no)) < 0.75*parameters.shock_val))
cell->set_coarsen_flag();
transfer_in.push_back(predictor);
triangulation.prepare_coarsening_and_refinement();
-
+
SolutionTransfer<dim, double> soltrans(dof_handler);
soltrans.prepare_for_coarsening_and_refinement(transfer_in);
template <int dim>
void ConservationLaw<dim>::output_results () const
{
- typename EulerEquations<dim>::Postprocessor
+ typename EulerEquations<dim>::Postprocessor
postprocessor (parameters.schlieren_plot);
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
-
+
data_out.add_data_vector (current_solution,
EulerEquations<dim>::component_names (),
DataOut<dim>::type_dof_data,
// solution. At the end of this process, we
// output the initial solution.
template <int dim>
-void ConservationLaw<dim>::run ()
+void ConservationLaw<dim>::run ()
{
{
GridIn<dim> grid_in;
std::ifstream input_file(parameters.mesh_filename.c_str());
Assert (input_file, ExcFileNotOpen(parameters.mesh_filename.c_str()));
- grid_in.read_ucd(input_file);
+ grid_in.read_ucd(input_file);
}
-
+
dof_handler.clear();
dof_handler.distribute_dofs (fe);
-
+
// Size all of the fields.
old_solution.reinit (dof_handler.n_dofs());
current_solution.reinit (dof_handler.n_dofs());
compute_refinement_indicators(refinement_indicators);
refine_grid(refinement_indicators);
-
+
setup_system();
VectorTools::interpolate(dof_handler,
<< dof_handler.n_dofs()
<< std::endl
<< std::endl;
-
+
std::cout << " NonLin Res Lin Iter Lin Res" << std::endl
<< " _____________________________________" << std::endl;
{
Matrix->PutScalar(0);
Matrix->FillComplete();
-
+
right_hand_side = 0;
assemble_system ();
-
+
const double res_norm = right_hand_side.l2_norm();
if (std::fabs(res_norm) < 1e-10)
{
std::pair<unsigned int, double> convergence
= solve (newton_update);
-
+
current_solution += newton_update;
-
+
std::printf(" %-16.3e %04d %-5.2e\n",
res_norm, convergence.first, convergence.second);
}
++nonlin_iter;
AssertThrow (nonlin_iter <= 10,
ExcMessage ("No convergence in nonlinear solver"));
- }
+ }
// We only get to this point if the
// Newton iteration has converged, so
// that the program aborts if no
// input file name is given on the
// command line.
-int main (int argc, char *argv[])
+int main (int argc, char *argv[])
{
if (argc != 2)
{
std::cout << "Usage:" << argv[0] << " infile" << std::endl;
std::exit(1);
}
-
+
try
{
ConservationLaw<2> cons (argv[1]);
<< std::endl;
return 1;
}
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
};
-
+
return 0;
}
projects / dealii-svn.git / commitdiff