The source file coupled_system.h with comments:
#ifndef __coupled_system_h__
#define __coupled_system_h__
DiffSystem framework files
#include "libmesh/fem_function_base.h"
#include "libmesh/fem_system.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/parameter_vector.h"
using namespace libMesh;
The Navier-Stokes system class.
FEMSystem, TimeSolver and NewtonSolver will handle most tasks,
but we must specify element residuals
class CoupledSystem : public FEMSystem
{
public:
Constructor
CoupledSystem(EquationSystems& es,
const std::string& name_in,
const unsigned int number_in)
: FEMSystem(es, name_in, number_in), Peclet(1.) {qoi.resize(1);}
Function to get computed QoI values
Number &get_QoI_value()
{
return computed_QoI;
}
Number &get_parameter_value(unsigned int parameter_index)
{
return parameters[parameter_index];
}
ParameterVector &get_parameter_vector()
{
parameter_vector.resize(parameters.size());
for(unsigned int i = 0; i != parameters.size(); ++i)
{
parameter_vector[i] = ¶meters[i];
}
return parameter_vector;
}
Real &get_Pe()
{
return Peclet;
}
protected:
System initialization
virtual void init_data ();
Context initialization
virtual void init_context(DiffContext &context);
Element residual and jacobian calculations
Time dependent parts
virtual bool element_time_derivative (bool request_jacobian,
DiffContext& context);
Constraint parts
virtual bool element_constraint (bool request_jacobian,
DiffContext& context);
Postprocessed output
virtual void postprocess ();
Parameters associated with the system
std::vector<Number> parameters;
Indices for each variable;
unsigned int p_var, u_var, v_var, C_var;
The ParameterVector object that will contain pointers to
the system parameters
ParameterVector parameter_vector;
The Peclet number for the species transport
Real Peclet;
The functionals to be computed as QoIs
Number computed_QoI;
};
class CoupledFEMFunctionsx : public FEMFunctionBase<Number>
{
public:
Constructor
CoupledFEMFunctionsx(System& /* sys */, unsigned int var_number) {var = var_number;}
Destructor
virtual ~CoupledFEMFunctionsx () {}
virtual AutoPtr<FEMFunctionBase<Number> > clone () const
{return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsx(*this) ); }
virtual void operator() (const FEMContext&, const Point&,
const Real, DenseVector<Number>&)
{libmesh_error();}
virtual Number operator() (const FEMContext&, const Point& p,
const Real time = 0.);
private:
unsigned int var;
};
class CoupledFEMFunctionsy : public FEMFunctionBase<Number>
{
public:
Constructor
CoupledFEMFunctionsy(System& /* sys */, unsigned int var_number) {var = var_number;}
Destructor
virtual ~CoupledFEMFunctionsy () {}
virtual AutoPtr<FEMFunctionBase<Number> > clone () const
{return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsy(*this) ); }
virtual void operator() (const FEMContext&, const Point&,
const Real,
DenseVector<Number>&)
{libmesh_error();}
virtual Number operator() (const FEMContext&, const Point& p,
const Real time = 0.);
private:
unsigned int var;
};
#endif //__coupled_system_h__
The source file domain.h with comments:
#include "femparameters.h"
#include "libmesh/mesh.h"
Bring in everything from the libMesh namespace
class FEMParameters;
void build_domain (Mesh &mesh, FEMParameters ¶m);
The source file femparameters.h with comments:
#ifndef __fem_parameters_h__
#define __fem_parameters_h__
#include <limits>
#include "libmesh/libmesh_common.h"
#include "libmesh/dof_map.h"
#include "libmesh/enum_norm_type.h"
#include "libmesh/getpot.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
class FEMParameters
{
public:
FEMParameters() :
initial_timestep(0), n_timesteps(100),
transient(true),
deltat_reductions(0),
timesolver_core("euler"),
end_time(std::numeric_limits<Real>::max()),
deltat(0.0001), timesolver_theta(0.5),
timesolver_maxgrowth(0.), timesolver_tolerance(0.),
timesolver_upper_tolerance(0.),
steadystate_tolerance(0.),
timesolver_norm(0,L2),
dimension(2),
domaintype("square"), domainfile("mesh.xda"), elementtype("quad"),
fine_mesh_file_primal("fine_mesh.xda"), fine_mesh_soln_primal("fine_mesh_soln.xda"),
fine_mesh_file_adjoint("fine_mesh.xda"), fine_mesh_soln_adjoint("fine_mesh_soln.xda"),
elementorder(2),
domain_xmin(0.0), domain_ymin(0.0), domain_zmin(0.0),
domain_edge_width(1.0), domain_edge_length(1.0), domain_edge_height(1.0),
coarsegridx(1), coarsegridy(1), coarsegridz(1),
coarserefinements(0), coarsecoarsenings(0), extrarefinements(0),
use_petsc_snes(false),
time_solver_quiet(true), solver_quiet(true),
reuse_preconditioner(false),
require_residual_reduction(true),
min_step_length(1e-5),
max_linear_iterations(200000), max_nonlinear_iterations(20),
relative_step_tolerance(1.e-7), relative_residual_tolerance(1.e-10),
initial_linear_tolerance(1.e-3), minimum_linear_tolerance(TOLERANCE*TOLERANCE),
linear_tolerance_multiplier(1.e-3),
nelem_target(30000), global_tolerance(0.0),
refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
refine_uniformly(false),
max_adaptivesteps(1),
initial_adaptivesteps(0), initial_sobolev_order(1),
initial_extra_quadrature(0),
indicator_type("kelly"), patch_reuse(true), sobolev_order(1), adjoint_residual_type("patchpatch"), alternate_with_uniform_steps("false"), alternate_step_number(2),
component_wise_error("false"), compare_to_fine_solution_primal(false), compare_to_fine_solution_adjoint(false), compute_sensitivities(false), do_forward_sensitivity(true), do_adjoint_sensitivity(false),
postprocess_adjoint(false),
write_interval(10),
write_gmv_error(false), write_tecplot_error(false),
output_xda(false), output_xdr(false),
output_bz2(true), output_gz(true),
output_gmv(false), output_tecplot(false),
run_simulation(true), run_postprocess(false),
fe_family(1, "LAGRANGE"), fe_order(1, 1),
extra_quadrature_order(0),
analytic_jacobians(true), verify_analytic_jacobians(0.0),
print_solution_norms(false), print_solutions(false),
print_residual_norms(false), print_residuals(false),
print_jacobian_norms(false), print_jacobians(false) {}
void read(GetPot &input);
unsigned int initial_timestep, n_timesteps;
bool transient;
unsigned int deltat_reductions;
std::string timesolver_core;
Real end_time, deltat, timesolver_theta,
timesolver_maxgrowth, timesolver_tolerance,
timesolver_upper_tolerance, steadystate_tolerance;
std::vector<FEMNormType> timesolver_norm;
unsigned int dimension;
std::string domaintype, domainfile, elementtype, fine_mesh_file_primal, fine_mesh_soln_primal, fine_mesh_file_adjoint, fine_mesh_soln_adjoint;
Real elementorder;
Real domain_xmin, domain_ymin, domain_zmin;
Real domain_edge_width, domain_edge_length, domain_edge_height;
unsigned int coarsegridx, coarsegridy, coarsegridz;
unsigned int coarserefinements, coarsecoarsenings, extrarefinements;
bool use_petsc_snes;
bool time_solver_quiet, solver_quiet, reuse_preconditioner, require_residual_reduction;
Real min_step_length;
unsigned int max_linear_iterations, max_nonlinear_iterations;
Real relative_step_tolerance, relative_residual_tolerance,
initial_linear_tolerance, minimum_linear_tolerance,
linear_tolerance_multiplier;
unsigned int nelem_target;
Real global_tolerance;
Real refine_fraction, coarsen_fraction, coarsen_threshold;
bool refine_uniformly;
unsigned int max_adaptivesteps;
unsigned int initial_adaptivesteps;
unsigned int initial_sobolev_order;
unsigned int initial_extra_quadrature;
std::string indicator_type;
bool patch_reuse;
unsigned int sobolev_order;
std::string adjoint_residual_type;
bool alternate_with_uniform_steps;
unsigned int alternate_step_number;
bool component_wise_error;
bool compare_to_fine_solution_primal, compare_to_fine_solution_adjoint;
bool compute_sensitivities;
bool do_forward_sensitivity;
bool do_adjoint_sensitivity;
bool postprocess_adjoint;
unsigned int write_interval;
bool write_gmv_error, write_tecplot_error, output_xda, output_xdr,
output_bz2, output_gz, output_gmv, output_tecplot;
bool run_simulation, run_postprocess;
std::vector<std::string> fe_family;
std::vector<unsigned int> fe_order;
int extra_quadrature_order;
bool analytic_jacobians;
Real verify_analytic_jacobians;
bool print_solution_norms, print_solutions,
print_residual_norms, print_residuals,
print_jacobian_norms, print_jacobians;
};
#endif // __fem_parameters_h__
The source file H-qoi.h with comments:
#ifndef H_QOI_H
#define H_QOI_H
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
#include "libmesh/diff_qoi.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
class CoupledSystemQoI : public DifferentiableQoI
{
public:
CoupledSystemQoI(){}
virtual ~CoupledSystemQoI(){}
virtual void init_qoi( std::vector<Number>& sys_qoi);
virtual void postprocess( ){}
virtual void side_qoi_derivative(DiffContext &context, const QoISet & qois);
virtual void side_qoi(DiffContext &context, const QoISet & qois);
virtual AutoPtr<DifferentiableQoI> clone( )
{ AutoPtr<DifferentiableQoI> my_clone( new CoupledSystemQoI );
*my_clone = *this;
return my_clone;
}
};
#endif // H_QOI_H
The source file initial.h with comments:
#include "libmesh/parameters.h"
#include "libmesh/point.h"
#include "libmesh/vector_value.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
void read_initial_parameters();
void finish_initialization();
Number initial_value(const Point& p,
const Parameters&,
const std::string&,
const std::string&);
Gradient initial_grad(const Point& p,
const Parameters&,
const std::string&,
const std::string&);
The source file mysystems.h with comments:
#include "coupled_system.h"
#include "femparameters.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
FEMSystem &build_system(EquationSystems &es, GetPot &, FEMParameters& /* param */)
{
CoupledSystem &system = es.add_system<CoupledSystem> ("CoupledSystem");
return system;
}
Number exact_value(const Point& /* p */, // xyz location
const Parameters& /* param */, // EquationSystem parameters
const std::string&, // sys_name
const std::string&) // unknown_name
{
std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
return 0;
}
Gradient exact_grad(const Point& /* p */, // xyz location
const Parameters& /* param */, // EquationSystems parameters
const std::string&, // sys_name
const std::string&) // unknown_name
{
std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
return 0;
}
The source file output.h with comments:
#include <string>
Bring in everything from the libMesh namespace
using namespace libMesh;
void start_output(unsigned int timesteps,
std::string filename,
std::string varname);
The source file adjoints_ex3.C with comments:
#include <iostream>
#include <sys/time.h>
#include <iomanip>
Libmesh includes
#include "libmesh/equation_systems.h"
#include "libmesh/twostep_time_solver.h"
#include "libmesh/euler_solver.h"
#include "libmesh/euler2_solver.h"
#include "libmesh/steady_solver.h"
#include "libmesh/newton_solver.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/petsc_diff_solver.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_tools.h"
#include "libmesh/mesh_base.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/system.h"
#include "libmesh/system_norm.h"
#include "libmesh/adjoint_residual_error_estimator.h"
#include "libmesh/const_fem_function.h"
#include "libmesh/error_vector.h"
#include "libmesh/fem_function_base.h"
#include "libmesh/getpot.h"
#include "libmesh/gmv_io.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/parameter_vector.h"
#include "libmesh/patch_recovery_error_estimator.h"
#include "libmesh/petsc_vector.h"
#include "libmesh/sensitivity_data.h"
#include "libmesh/tecplot_io.h"
#include "libmesh/uniform_refinement_estimator.h"
#include "libmesh/qoi_set.h"
#include "libmesh/weighted_patch_recovery_error_estimator.h"
Local includes
#include "coupled_system.h"
#include "domain.h"
#include "initial.h"
#include "femparameters.h"
#include "mysystems.h"
#include "output.h"
#include "H-qoi.h"
We solve a coupled Stokes + Convection Diffusion system in an H channel geometry
with 2 inlets and 2 outlets. The QoI is the species
flux from left outlet
Channel Geometry: Wall ---------------------------------------------------------------- 0 1 ---------------------------- ------------------------------- | | Wall | | Wall | | ----------------------------- ------------------------------- 2 2 ----------------------------------------------------------------- Wall
The equations governing this flow are: Stokes: -VectorLaplacian(velocity) + grad(pressure) = vector(0) Convection-Diffusion: - dot(velocity, grad(concentration) ) + Laplacian(concentration) = 0
The boundary conditions are: u_1(0) = -(y-2)*(y-3), u_2(0) = 0 ; u_1(1) = (y-2)*(y-3), u_2(1) = 0 ; u_1(walls) = 0, u_2(walls) = 0; C(0) = 1 ; C(1) = 0; grad(C) dot n (walls) = 0 ; grad(C) dot n (2) = 0 ; grad(C) dot n (3) = 0
The QoI is: Q((u,v), C) = integral_{left_outlet} - u * C ds
The complete equal order adjoint QoI error estimate is: (Notation for derivatives: grad(C) = C,1 + C,2) Q(u) - Q(u_h) \leq |e(u_1)|_{H1} |e(u_1^*)|_{H1} + |e(u_2)|_{H1} |e(u_2^*)|_{H1} + (1/Pe) * |e(C)|_{H1} |e(C^*)|_{H1} + ||e(u_1,1)||_{L2} ||e(p^*)||_{L2} + ||e(u_2,2)||_{L2} ||e(p^*)||_{L2} + ||e(u_1,1^*)||_{L2} ||e(p)||_{L2} + ||e(u_2,2^*)||_{L2} ||e(p)||_{L2} + ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} + ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2} ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} + ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2} = error_non_pressure + error_with_pressure + error_convection_diffusion_x + error_convection_diffusion_y
Bring in everything from the libMesh namespace
Channel Geometry: Wall ---------------------------------------------------------------- 0 1 ---------------------------- ------------------------------- | | Wall | | Wall | | ----------------------------- ------------------------------- 2 2 ----------------------------------------------------------------- Wall
The equations governing this flow are: Stokes: -VectorLaplacian(velocity) + grad(pressure) = vector(0) Convection-Diffusion: - dot(velocity, grad(concentration) ) + Laplacian(concentration) = 0
The boundary conditions are: u_1(0) = -(y-2)*(y-3), u_2(0) = 0 ; u_1(1) = (y-2)*(y-3), u_2(1) = 0 ; u_1(walls) = 0, u_2(walls) = 0; C(0) = 1 ; C(1) = 0; grad(C) dot n (walls) = 0 ; grad(C) dot n (2) = 0 ; grad(C) dot n (3) = 0
The QoI is: Q((u,v), C) = integral_{left_outlet} - u * C ds
The complete equal order adjoint QoI error estimate is: (Notation for derivatives: grad(C) = C,1 + C,2) Q(u) - Q(u_h) \leq |e(u_1)|_{H1} |e(u_1^*)|_{H1} + |e(u_2)|_{H1} |e(u_2^*)|_{H1} + (1/Pe) * |e(C)|_{H1} |e(C^*)|_{H1} + ||e(u_1,1)||_{L2} ||e(p^*)||_{L2} + ||e(u_2,2)||_{L2} ||e(p^*)||_{L2} + ||e(u_1,1^*)||_{L2} ||e(p)||_{L2} + ||e(u_2,2^*)||_{L2} ||e(p)||_{L2} + ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} + ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2} ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} + ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2} = error_non_pressure + error_with_pressure + error_convection_diffusion_x + error_convection_diffusion_y
Bring in everything from the libMesh namespace
using namespace libMesh;
Number output files
std::string numbered_filename(unsigned int t_step, // The timestep count
unsigned int a_step, // The adaptive step count
std::string solution_type, // primal or adjoint solve
std::string type,
std::string extension,
FEMParameters ¶m)
{
std::ostringstream file_name;
file_name << solution_type
<< ".out."
<< type
<< '.'
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step
<< '.'
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step
<< '.'
<< extension;
if (param.output_bz2)
file_name << ".bz2";
else if (param.output_gz)
file_name << ".gz";
return file_name.str();
}
Write tecplot, gmv and xda/xdr output
void write_output(EquationSystems &es,
unsigned int t_step, // The timestep count
unsigned int a_step, // The adaptive step count
std::string solution_type, // primal or adjoint solve
FEMParameters ¶m)
{
MeshBase &mesh = es.get_mesh();
#ifdef LIBMESH_HAVE_GMV
if (param.output_gmv)
{
std::ostringstream file_name_gmv;
file_name_gmv << solution_type
<< ".out.gmv."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step
<< '.'
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step;
GMVIO(mesh).write_equation_systems
(file_name_gmv.str(), es);
}
#endif
#ifdef LIBMESH_HAVE_TECPLOT_API
if (param.output_tecplot)
{
std::ostringstream file_name_tecplot;
file_name_tecplot << solution_type
<< ".out."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step
<< '.'
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step
<< ".plt";
TecplotIO(mesh).write_equation_systems
(file_name_tecplot.str(), es);
}
#endif
if (param.output_xda || param.output_xdr)
mesh.renumber_nodes_and_elements();
if (param.output_xda)
{
mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param));
es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xda", param),
libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
EquationSystems::WRITE_ADDITIONAL_DATA);
}
if (param.output_xdr)
{
mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param));
es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param),
libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
EquationSystems::WRITE_ADDITIONAL_DATA);
}
}
void write_output_headers(FEMParameters ¶m)
{
Only one processor needs to take care of headers.
if (libMesh::processor_id() != 0)
return;
start_output(param.initial_timestep, "out_clocktime.m", "vector_clocktime");
if (param.run_simulation)
{
start_output(param.initial_timestep, "out_activemesh.m", "table_activemesh");
start_output(param.initial_timestep, "out_solvesteps.m", "table_solvesteps");
if (param.timesolver_tolerance)
{
start_output(param.initial_timestep, "out_time.m", "vector_time");
start_output(param.initial_timestep, "out_timesteps.m", "vector_timesteps");
}
if (param.steadystate_tolerance)
start_output(param.initial_timestep, "out_changerate.m", "vector_changerate");
}
}
void write_output_solvedata(EquationSystems &es,
unsigned int a_step,
unsigned int newton_steps,
unsigned int krylov_steps,
unsigned int tv_sec,
unsigned int tv_usec)
{
MeshBase &mesh = es.get_mesh();
unsigned int n_active_elem = mesh.n_active_elem();
unsigned int n_active_dofs = es.n_active_dofs();
if (libMesh::processor_id() == 0)
{
Write out the number of elements/dofs used
std::ofstream activemesh ("out_activemesh.m",
std::ios_base::app | std::ios_base::out);
activemesh.precision(17);
activemesh << (a_step + 1) << ' '
<< n_active_elem << ' '
<< n_active_dofs << std::endl;
Write out the number of solver steps used
std::ofstream solvesteps ("out_solvesteps.m",
std::ios_base::app | std::ios_base::out);
solvesteps.precision(17);
solvesteps << newton_steps << ' '
<< krylov_steps << std::endl;
Write out the clock time
std::ofstream clocktime ("out_clocktime.m",
std::ios_base::app | std::ios_base::out);
clocktime.precision(17);
clocktime << tv_sec << '.' << tv_usec << std::endl;
}
}
void write_output_footers(FEMParameters ¶m)
{
Write footers on output .m files
if (libMesh::processor_id() == 0)
{
std::ofstream clocktime ("out_clocktime.m",
std::ios_base::app | std::ios_base::out);
clocktime << "];" << std::endl;
if (param.run_simulation)
{
std::ofstream activemesh ("out_activemesh.m",
std::ios_base::app | std::ios_base::out);
activemesh << "];" << std::endl;
std::ofstream solvesteps ("out_solvesteps.m",
std::ios_base::app | std::ios_base::out);
solvesteps << "];" << std::endl;
if (param.timesolver_tolerance)
{
std::ofstream times ("out_time.m",
std::ios_base::app | std::ios_base::out);
times << "];" << std::endl;
std::ofstream timesteps ("out_timesteps.m",
std::ios_base::app | std::ios_base::out);
timesteps << "];" << std::endl;
}
if (param.steadystate_tolerance)
{
std::ofstream changerate ("out_changerate.m",
std::ios_base::app | std::ios_base::out);
changerate << "];" << std::endl;
}
}
We're done, so write out a file (for e.g. Make to check)
std::ofstream complete ("complete");
complete << "complete" << std::endl;
}
}
#if defined(LIBMESH_HAVE_GMV) || defined(LIBMESH_HAVE_TECPLOT_API)
void write_error(EquationSystems &es,
ErrorVector &error,
unsigned int t_number,
unsigned int a_number,
FEMParameters ¶m,
std::string error_type)
#else
void write_error(EquationSystems &,
ErrorVector &,
unsigned int,
unsigned int,
FEMParameters &,
std::string)
#endif
{
#ifdef LIBMESH_HAVE_GMV
if (param.write_gmv_error)
{
std::ostringstream error_gmv;
error_gmv << "error.gmv."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< a_number
<< "."
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< t_number
<< error_type;
error.plot_error(error_gmv.str(), es.get_mesh());
}
#endif
#ifdef LIBMESH_HAVE_TECPLOT_API
if (param.write_tecplot_error)
{
std::ostringstream error_tecplot;
error_tecplot << "error.plt."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< a_number
<< "."
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< t_number
<< error_type;
error.plot_error(error_tecplot.str(), es.get_mesh());
}
#endif
}
void read_output(EquationSystems &es,
unsigned int t_step,
unsigned int a_step,
std::string solution_type,
FEMParameters ¶m)
{
MeshBase &mesh = es.get_mesh();
std::string file_name_mesh, file_name_soln;
Look for ASCII files first
if (param.output_xda)
{
file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param);
file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xda", param);
}
else if (param.output_xdr)
{
file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param);
file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param);
}
Read in the mesh
mesh.read(file_name_mesh);
And the stored solution
es.read(file_name_soln, libMeshEnums::READ,
EquationSystems::READ_HEADER |
EquationSystems::READ_DATA |
EquationSystems::READ_ADDITIONAL_DATA);
Put systems in a consistent state
for (unsigned int i = 0; i != es.n_systems(); ++i)
es.get_system<FEMSystem>(i).update();
Figure out the current time
Real current_time = 0., current_timestep = 0.;
if (param.timesolver_tolerance)
{
std::ifstream times ("out_time.m");
std::ifstream timesteps ("out_timesteps.m");
if (times.is_open() && timesteps.is_open())
{
Read headers
const unsigned int headersize = 25;
char header[headersize];
timesteps.getline (header, headersize);
if (strcmp(header, "vector_timesteps = [") != 0)
{
std::cout << "Bad header in out_timesteps.m:" << std::endl
<< header
<< std::endl;
libmesh_error();
}
times.getline (header, headersize);
if (strcmp(header, "vector_time = [") != 0)
{
std::cout << "Bad header in out_time.m:" << std::endl
<< header
<< std::endl;
libmesh_error();
}
Read each timestep
for (unsigned int i = 0; i != t_step; ++i)
{
if (!times.good())
libmesh_error();
times >> current_time;
timesteps >> current_timestep;
}
Remember to increment the last timestep; out_times.m
lists each *start* time
current_time += current_timestep;
}
else
libmesh_error();
}
else
current_time = t_step * param.deltat;
for (unsigned int i = 0; i != es.n_systems(); ++i)
es.get_system<FEMSystem>(i).time = current_time;
}
void set_system_parameters(FEMSystem &system, FEMParameters ¶m)
{
Verify analytic jacobians against numerical ones?
system.verify_analytic_jacobians = param.verify_analytic_jacobians;
More desperate debugging options
system.print_solution_norms = param.print_solution_norms;
system.print_solutions = param.print_solutions;
system.print_residual_norms = param.print_residual_norms;
system.print_residuals = param.print_residuals;
system.print_jacobian_norms = param.print_jacobian_norms;
system.print_jacobians = param.print_jacobians;
Solve this as a time-dependent or steady system
if (param.transient)
{
UnsteadySolver *innersolver;
if (param.timesolver_core == "euler2")
{
Euler2Solver *euler2solver =
new Euler2Solver(system);
euler2solver->theta = param.timesolver_theta;
innersolver = euler2solver;
}
else if (param.timesolver_core == "euler")
{
EulerSolver *eulersolver =
new EulerSolver(system);
eulersolver->theta = param.timesolver_theta;
innersolver = eulersolver;
}
else
{
std::cerr << "Don't recognize core TimeSolver type: "
<< param.timesolver_core << std::endl;
libmesh_error();
}
if (param.timesolver_tolerance)
{
TwostepTimeSolver *timesolver =
new TwostepTimeSolver(system);
timesolver->max_growth = param.timesolver_maxgrowth;
timesolver->target_tolerance = param.timesolver_tolerance;
timesolver->upper_tolerance = param.timesolver_upper_tolerance;
timesolver->component_norm = SystemNorm(param.timesolver_norm);
timesolver->core_time_solver =
AutoPtr<UnsteadySolver>(innersolver);
system.time_solver =
AutoPtr<UnsteadySolver>(timesolver);
}
else
system.time_solver =
AutoPtr<TimeSolver>(innersolver);
}
else
system.time_solver =
AutoPtr<TimeSolver>(new SteadySolver(system));
system.time_solver->reduce_deltat_on_diffsolver_failure =
param.deltat_reductions;
system.time_solver->quiet = param.time_solver_quiet;
Set the time stepping options
system.deltat = param.deltat;
And the integration options
system.extra_quadrature_order = param.extra_quadrature_order;
And the nonlinear solver options
if (param.use_petsc_snes)
{
#ifdef LIBMESH_HAVE_PETSC
PetscDiffSolver *solver = new PetscDiffSolver(system);
system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
#else
libmesh_error();
#endif
}
else
{
NewtonSolver *solver = new NewtonSolver(system);
system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
solver->quiet = param.solver_quiet;
solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
solver->minsteplength = param.min_step_length;
solver->relative_step_tolerance = param.relative_step_tolerance;
solver->relative_residual_tolerance = param.relative_residual_tolerance;
solver->require_residual_reduction = param.require_residual_reduction;
solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
if (system.time_solver->reduce_deltat_on_diffsolver_failure)
{
solver->continue_after_max_iterations = true;
solver->continue_after_backtrack_failure = true;
}
And the linear solver options
solver->max_linear_iterations = param.max_linear_iterations;
solver->initial_linear_tolerance = param.initial_linear_tolerance;
solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
}
}
#ifdef LIBMESH_ENABLE_AMR
AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
FEMParameters ¶m)
{
AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
mesh_refinement->coarsen_by_parents() = true;
mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
mesh_refinement->nelem_target() = param.nelem_target;
mesh_refinement->refine_fraction() = param.refine_fraction;
mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
return mesh_refinement;
}
#endif // LIBMESH_ENABLE_AMR
This function builds the Kelly error indicator. This indicator can be used
for comparisons of adjoint and non-adjoint based error indicators
AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters& /* param */)
{
AutoPtr<ErrorEstimator> error_estimator;
error_estimator.reset(new KellyErrorEstimator);
return error_estimator;
}
Functions to build the adjoint based error indicators
The error_non_pressure and error_pressure constributions are estimated using
the build_error_estimator_component_wise function below
AutoPtr<ErrorEstimator>
build_error_estimator_component_wise
(FEMParameters ¶m,
std::vector<std::vector<Real> > &term_weights,
std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type)
{
AutoPtr<ErrorEstimator> error_estimator;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
Both the primal and dual weights are going to be estimated using the patch recovery error estimator
PatchRecoveryErrorEstimator *p1 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->primal_error_estimator().reset(p1);
PatchRecoveryErrorEstimator *p2 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->dual_error_estimator().reset(p2);
Set the boolean for specifying whether we are reusing patches while building the patch recovery estimates
p1->set_patch_reuse(param.patch_reuse);
p2->set_patch_reuse(param.patch_reuse);
Using the user filled error norm type vector, we pass the type of norm to be used for
the error in each variable, we can have different types of norms for the primal and
dual variables
unsigned int size = primal_error_norm_type.size();
libmesh_assert_equal_to (size, dual_error_norm_type.size());
for (unsigned int i = 0; i != size; ++i)
{
adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
}
Now we set the right weights for each term in the error estimate, using the user provided
term_weights matrix
libmesh_assert_equal_to (size, term_weights.size());
for (unsigned int i = 0; i != term_weights.size(); ++i)
{
libmesh_assert_equal_to (size, term_weights[i].size());
adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
for (unsigned int j = 0; j != size; ++j)
if (i != j)
adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
}
return error_estimator;
}
The error_convection_diffusion_x and error_convection_diffusion_y are the nonlinear contributions which
are computed using the build_weighted_error_estimator_component_wise below
AutoPtr<ErrorEstimator>
build_weighted_error_estimator_component_wise
(FEMParameters ¶m,
std::vector<std::vector<Real> > &term_weights,
std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type,
std::vector<FEMFunctionBase<Number>*> coupled_system_weight_functions)
{
AutoPtr<ErrorEstimator> error_estimator;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
Using the user filled error norm type vector, we pass the type of norm to be used for
the error in each variable, we can have different types of norms for the primal and
dual variables
WeightedPatchRecoveryErrorEstimator *p1 =
new WeightedPatchRecoveryErrorEstimator;
adjoint_residual_estimator->primal_error_estimator().reset(p1);
PatchRecoveryErrorEstimator *p2 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->dual_error_estimator().reset(p2);
p1->set_patch_reuse(param.patch_reuse);
p2->set_patch_reuse(param.patch_reuse);
This is the critical difference with the build_error_estimate_component_wise
p1->weight_functions.clear();
We pass the pointers to the user specified weight functions to the patch recovery
error estimator objects declared above
unsigned int size = primal_error_norm_type.size();
libmesh_assert(coupled_system_weight_functions.size() == size);
libmesh_assert(dual_error_norm_type.size() == size);
for (unsigned int i = 0; i != size; ++i)
{
p1->weight_functions.push_back(coupled_system_weight_functions[i]);
adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
}
Now we set the right weights for each term in the error estimate, using the user provided
term_weights matrix
libmesh_assert_equal_to (size, term_weights.size());
for (unsigned int i = 0; i != term_weights.size(); ++i)
{
libmesh_assert_equal_to (size, term_weights[i].size());
adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
for (unsigned int j = 0; j != size; ++j)
if (i != j)
adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
}
return error_estimator;
}
The main program.
int main (int argc, char** argv)
{
Initialize libMesh.
LibMeshInit init (argc, argv);
Skip adaptive examples on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
std::cout << "Started " << argv[0] << std::endl;
Make sure the general input file exists, and parse it
{
std::ifstream i("general.in");
if (!i)
{
std::cerr << '[' << libMesh::processor_id()
<< "] Can't find general.in; exiting early."
<< std::endl;
libmesh_error();
}
}
GetPot infile("general.in");
Read in parameters from the input file
FEMParameters param;
param.read(infile);
Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Create a mesh.
Mesh mesh (param.dimension);
And an object to refine it
AutoPtr<MeshRefinement> mesh_refinement =
build_mesh_refinement(mesh, param);
And an EquationSystems to run on it
EquationSystems equation_systems (mesh);
std::cout << "Building mesh" << std::endl;
if (!param.initial_timestep && param.run_simulation)
build_domain(mesh, param);
std::cout << "Building system" << std::endl;
FEMSystem &system = build_system(equation_systems, infile, param);
set_system_parameters(system, param);
std::cout << "Initializing systems" << std::endl;
Create a QoI object, we will need this to compute the QoI
{
CoupledSystemQoI qoi;
Our QoI is computed on the side
qoi.assemble_qoi_sides = true;
system.attach_qoi(&qoi);
}
equation_systems.init ();
Print information about the mesh and system to the screen.
mesh.print_info();
equation_systems.print_info();
std::cout<<"Starting adaptive loop"<<std::endl<<std::endl;
Count the number of steps used
unsigned int a_step = 0;
Get the linear solver object to set the preconditioner reuse flag
LinearSolver<Number> *linear_solver = system.get_linear_solver();
Adaptively solve the timestep
for (; a_step <= param.max_adaptivesteps; ++a_step)
{
We have to ensure that we are refining to either an error
tolerance or to a target number of elements, and, not to a
non-zero tolerance and a target number of elements
simultaneously
if(param.global_tolerance > 0 && param.nelem_target > 0)
{
std::cout<<"We cant refine to both a non-zero tolerance and a target number of elements, EXITING adaptive loop. "<<std::endl<<std::endl;
break;
}
We dont want to reuse the preconditioner for the forward solve
linear_solver->reuse_preconditioner(false);
std::cout<< "Adaptive step " << a_step << std::endl;
std::cout<< "Solving the forward problem" <<std::endl;
std::cout << "We have " << mesh.n_active_elem()
<< " active elements and " << equation_systems.n_active_dofs()
<< " active dofs." << std::endl << std::endl;
Solve the forward system
system.solve();
Write the output
write_output(equation_systems, 0, a_step, "primal", param);
Compute the QoI
system.assemble_qoi();
We just call a postprocess here to set the variable computed_QoI to the value computed by assemble_qoi
system.postprocess();
Get the value of the computed_QoI variable of the CoupledSystem class
Number QoI_0_computed = (dynamic_cast<CoupledSystem&>(system)).get_QoI_value();
std::cout<< "The boundary QoI is " << std::setprecision(17) << QoI_0_computed << std::endl << std::endl;
You dont need to compute error estimates and refine at the last
adaptive step, only before that
if(a_step!=param.max_adaptivesteps)
{
NumericVector<Number> &primal_solution = (*system.solution);
Make sure we get the contributions to the adjoint RHS from the sides
system.assemble_qoi_sides = true;
We are about to solve the adjoint system, but before we
do this we see the same preconditioner flag to reuse the
preconditioner from the forward solver
linear_solver->reuse_preconditioner(param.reuse_preconditioner);
Solve the adjoint system
std::cout<< "Solving the adjoint problem" <<std::endl;
system.adjoint_solve();
Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unneccesarily in the error estimator
system.set_adjoint_already_solved(true);
To plot the adjoint solution, we swap it with the primal solution
and use the write_output function
NumericVector<Number> &dual_solution = system.get_adjoint_solution();
primal_solution.swap(dual_solution);
write_output(equation_systems, 0, a_step, "adjoint", param);
Swap back
primal_solution.swap(dual_solution);
We need the values of the parameters Pe from the
system for the adjoint error estimate
Real Pe = (dynamic_cast<CoupledSystem&>(system)).get_Pe();
The total error is the sum: error = error_non_pressure +
error_with_pressure + ...
error_estimator_convection_diffusion_x +
error_estimator_convection_diffusion_y
ErrorVector error;
We first construct the non-pressure contributions
ErrorVector error_non_pressure;
First we build the norm_type_vector_non_pressure vectors and
weights_matrix_non_pressure matrix for the non-pressure term
error contributions
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_non_pressure;
primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
primal_norm_type_vector_non_pressure.push_back(L2);
primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_non_pressure;
dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
dual_norm_type_vector_non_pressure.push_back(L2);
dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
std::vector<std::vector<Real> >
weights_matrix_non_pressure(system.n_vars(),
std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_non_pressure[0][0] = 1.;
weights_matrix_non_pressure[1][1] = 1.;
weights_matrix_non_pressure[3][3] = 1./Pe;
We build the error estimator to estimate the contributions
to the QoI error from the non pressure term
AutoPtr<ErrorEstimator> error_estimator_non_pressure =
build_error_estimator_component_wise
(param, weights_matrix_non_pressure,
primal_norm_type_vector_non_pressure,
dual_norm_type_vector_non_pressure);
Estimate the contributions to the QoI error from the non
pressure terms
error_estimator_non_pressure->estimate_error(system, error_non_pressure);
Plot the estimated error from the non_pressure terms
write_error(equation_systems, error_non_pressure, 0, a_step, param, "_non_pressure");
Now for the pressure contributions
ErrorVector error_with_pressure;
Next we build the norm_type_vector_with_pressure vectors and
weights_matrix_with_pressure matrix for the pressure term
error contributions
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_with_pressure;
primal_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
primal_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
primal_norm_type_vector_with_pressure.push_back(L2);
primal_norm_type_vector_with_pressure.push_back(L2);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_with_pressure;
dual_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
dual_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
dual_norm_type_vector_with_pressure.push_back(L2);
dual_norm_type_vector_with_pressure.push_back(L2);
std::vector<std::vector<Real> >
weights_matrix_with_pressure
(system.n_vars(),
std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_with_pressure[0][2] = 1.;
weights_matrix_with_pressure[1][2] = 1.;
weights_matrix_with_pressure[2][0] = 1.;
weights_matrix_with_pressure[2][1] = 1.;
We build the error estimator to estimate the contributions
to the QoI error from the pressure term
AutoPtr<ErrorEstimator> error_estimator_with_pressure =
build_error_estimator_component_wise
(param, weights_matrix_with_pressure,
primal_norm_type_vector_with_pressure,
dual_norm_type_vector_with_pressure);
Estimate the contributions to the QoI error from the pressure terms
error_estimator_with_pressure->estimate_error(system, error_with_pressure);
Plot the error due to the pressure terms
write_error(equation_systems, error_with_pressure, 0, a_step, param, "_with_pressure");
Now for the convection diffusion term errors (in the x and y directions)
ErrorVector error_convection_diffusion_x;
The norm type vectors and weights matrix for the convection_diffusion_x errors
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_convection_diffusion_x;
primal_norm_type_vector_convection_diffusion_x.push_back(L2);
primal_norm_type_vector_convection_diffusion_x.push_back(L2);
primal_norm_type_vector_convection_diffusion_x.push_back(L2);
primal_norm_type_vector_convection_diffusion_x.push_back(H1_X_SEMINORM);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_convection_diffusion_x;
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
Note that we need the error of the dual concentration in L2
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
std::vector<std::vector<Real> >
weights_matrix_convection_diffusion_x
(system.n_vars(),
std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_convection_diffusion_x[0][3] = 1.;
weights_matrix_convection_diffusion_x[3][3] = 1.;
We will also have to build and pass the weight functions to the weighted patch recovery estimators
We pass the identity function as weights to error entries that the above matrix will scale to 0.
We pass the identity function as weights to error entries that the above matrix will scale to 0.
ConstFEMFunction<Number> identity(1);
Declare object of class CoupledFEMFunctionsx, the definition of the function contains the weight
to be applied to the relevant terms
For ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2} term, returns C,1_h
CoupledFEMFunctionsx convdiffx0(system, 0);
For ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} term, returns (u_1)_h
CoupledFEMFunctionsx convdiffx3(system, 3);
Make a vector of pointers to these objects
std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_x;
coupled_system_weight_functions_x.push_back(&convdiffx0);
coupled_system_weight_functions_x.push_back(&identity);
coupled_system_weight_functions_x.push_back(&identity);
coupled_system_weight_functions_x.push_back(&convdiffx3);
Build the error estimator to estimate the contributions
to the QoI error from the convection diffusion x term
AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_x =
build_weighted_error_estimator_component_wise
(param, weights_matrix_convection_diffusion_x,
primal_norm_type_vector_convection_diffusion_x,
dual_norm_type_vector_convection_diffusion_x,
coupled_system_weight_functions_x);
Estimate the contributions to the QoI error from the
convection diffusion x term
error_estimator_convection_diffusion_x->estimate_error
(system, error_convection_diffusion_x);
Plot this error
write_error(equation_systems, error_convection_diffusion_x,
0, a_step, param, "_convection_diffusion_x");
Now for the y direction terms
ErrorVector error_convection_diffusion_y;
The norm type vectors and weights matrix for the convection_diffusion_x errors
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_convection_diffusion_y;
primal_norm_type_vector_convection_diffusion_y.push_back(L2);
primal_norm_type_vector_convection_diffusion_y.push_back(L2);
primal_norm_type_vector_convection_diffusion_y.push_back(L2);
primal_norm_type_vector_convection_diffusion_y.push_back(H1_Y_SEMINORM);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_convection_diffusion_y;
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
Note that we need the error of the dual concentration in L2
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
std::vector<std::vector<Real> >
weights_matrix_convection_diffusion_y
(system.n_vars(), std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_convection_diffusion_y[1][3] = 1.;
weights_matrix_convection_diffusion_y[3][3] = 1.;
For ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2} term, returns C,2_h
CoupledFEMFunctionsy convdiffy1(system, 1);
For ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} term, returns (u_2)_h
CoupledFEMFunctionsy convdiffy3(system, 3);
Make a vector of pointers to these objects
std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_y;
coupled_system_weight_functions_y.push_back(&identity);
coupled_system_weight_functions_y.push_back(&convdiffy1);
coupled_system_weight_functions_y.push_back(&identity);
coupled_system_weight_functions_y.push_back(&convdiffy3);
Build the error estimator to estimate the contributions
to the QoI error from the convection diffsion y term
AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_y =
build_weighted_error_estimator_component_wise
(param, weights_matrix_convection_diffusion_y,
primal_norm_type_vector_convection_diffusion_y,
dual_norm_type_vector_convection_diffusion_y,
coupled_system_weight_functions_y);
Estimate the contributions to the QoI error from the
convection diffusion y terms
error_estimator_convection_diffusion_y->estimate_error(system, error_convection_diffusion_y);
Plot this error
write_error(equation_systems, error_convection_diffusion_y, 0, a_step, param, "_convection_diffusion_y");
if(param.indicator_type == "adjoint_residual")
{
error.resize(error_non_pressure.size());
Now combine the contribs from the pressure and non
pressure terms to get the complete estimate
for(unsigned int i = 0; i < error.size(); i++)
{
error[i] = error_non_pressure[i] + error_with_pressure[i] + error_convection_diffusion_x[i] + error_convection_diffusion_y[i];
}
}
else
{
std::cout<<"Using Kelly Estimator"<<std::endl;
Build the Kelly error estimator
AutoPtr<ErrorEstimator> error_estimator = build_error_estimator(param);
Estimate the error
error_estimator->estimate_error(system, error);
}
write_error(equation_systems, error, 0, a_step, param, "_total");
We have to refine either based on reaching an error
tolerance or a number of elements target, which should be
verified above
Otherwise we flag elements by error tolerance or nelem target
Otherwise we flag elements by error tolerance or nelem target
if(param.refine_uniformly)
{
mesh_refinement->uniformly_refine(1);
}
else if(param.global_tolerance >= 0. && param.nelem_target == 0.) // Refine based on reaching an error tolerance
{
mesh_refinement->flag_elements_by_error_tolerance (error);
Carry out the adaptive mesh refinement/coarsening
mesh_refinement->refine_and_coarsen_elements();
}
else // Refine based on reaching a target number of elements
{
If we have reached the desired number of elements, we
dont do any more adaptive steps
if (mesh.n_active_elem() >= param.nelem_target)
{
std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
break;
}
mesh_refinement->flag_elements_by_nelem_target (error);
Carry out the adaptive mesh refinement/coarsening
mesh_refinement->refine_and_coarsen_elements();
}
equation_systems.reinit();
}
}
write_output_footers(param);
#endif // #ifndef LIBMESH_ENABLE_AMR
All done.
return 0;
}
The source file coupled_system.C with comments:
#include "libmesh/boundary_info.h"
#include "libmesh/dirichlet_boundaries.h"
#include "libmesh/dof_map.h"
#include "libmesh/fe_base.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fem_context.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh.h"
#include "libmesh/parallel.h"
#include "libmesh/quadrature.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/zero_function.h"
#include "coupled_system.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Function to set the Dirichlet boundary function for the inlet boundary velocities
class BdyFunction : public FunctionBase<Number>
{
public:
BdyFunction (unsigned int u_var, unsigned int v_var, int sign)
: _u_var(u_var), _v_var(v_var), _sign(sign)
{ this->_initialized = true; }
virtual Number operator() (const Point&, const Real = 0)
{ libmesh_not_implemented(); }
virtual void operator() (const Point& p,
const Real,
DenseVector<Number>& output)
{
output.resize(2);
output.zero();
const Real y=p(1);
Set the parabolic inflow boundary conditions at stations 0 & 1
output(_u_var) = (_sign)*((y-2) * (y-3));
output(_v_var) = 0;
}
virtual AutoPtr<FunctionBase<Number> > clone() const
{ return AutoPtr<FunctionBase<Number> > (new BdyFunction(_u_var, _v_var, _sign)); }
private:
const unsigned int _u_var, _v_var;
const Real _sign;
};
void CoupledSystem::init_data ()
{
Check the input file for Reynolds number, application type,
approximation type
GetPot infile("coupled_system.in");
Peclet = infile("Peclet", 1.);
unsigned int pressure_p = infile("pressure_p", 1);
std::string fe_family = infile("fe_family", std::string("LAGRANGE"));
LBB needs better-than-quadratic velocities for better-than-linear
pressures, and libMesh needs non-Lagrange elements for
better-than-quadratic velocities.
libmesh_assert((pressure_p == 1) || (fe_family != "LAGRANGE"));
FEFamily fefamily = Utility::string_to_enum<FEFamily>(fe_family);
Add the velocity components "u" & "v". They
will be approximated using second-order approximation.
u_var = this->add_variable ("u", static_cast<Order>(pressure_p+1),
fefamily);
v_var = this->add_variable ("v", static_cast<Order>(pressure_p+1),
fefamily);
Add the pressure variable "p". This will
be approximated with a first-order basis,
providing an LBB-stable pressure-velocity pair.
p_var = this->add_variable ("p", static_cast<Order>(pressure_p),
fefamily);
Add the Concentration variable "C". They will
be approximated using second-order approximation, the same as the velocity components
C_var = this->add_variable ("C", static_cast<Order>(pressure_p+1),
fefamily);
Tell the system to march velocity forward in time, but
leave p as a constraint only
this->time_evolving(u_var);
this->time_evolving(v_var);
this->time_evolving(C_var);
Useful debugging options
this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
this->print_jacobians = infile("print_jacobians", false);
this->print_element_jacobians = infile("print_element_jacobians", false);
Set Dirichlet boundary conditions
const boundary_id_type left_inlet_id = 0;
std::set<boundary_id_type> left_inlet_bdy;
left_inlet_bdy.insert(left_inlet_id);
const boundary_id_type right_inlet_id = 1;
std::set<boundary_id_type> right_inlet_bdy;
right_inlet_bdy.insert(right_inlet_id);
const boundary_id_type outlets_id = 2;
std::set<boundary_id_type> outlets_bdy;
outlets_bdy.insert(outlets_id);
const boundary_id_type wall_id = 3;
std::set<boundary_id_type> wall_bdy;
wall_bdy.insert(wall_id);
The uv identifier for the setting the inlet and wall velocity boundary conditions
std::vector<unsigned int> uv(1, u_var);
uv.push_back(v_var);
The C_only identifier for setting the concentrations at the inlets
std::vector<unsigned int> C_only(1, C_var);
The zero and constant functions
ZeroFunction<Number> zero;
ConstFunction<Number> one(1);
We need two boundary functions for the inlets, because the signs on the velocities
will be different
int velocity_sign = 1;
BdyFunction inflow_left(u_var, v_var, -velocity_sign);
BdyFunction inflow_right(u_var, v_var, velocity_sign);
On the walls we will apply the no slip and no penetration boundary condition, u=0, v=0
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (wall_bdy, uv, &zero));
On the inlet (left), we apply parabolic inflow boundary conditions for the velocity, u = - (y-2)*(y-3), v=0
and set C = 1
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (left_inlet_bdy, uv, &inflow_left));
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (left_inlet_bdy, C_only, &one));
On the inlet (right), we apply parabolic inflow boundary conditions for the velocity, u = (y-2)*(y-3), v=0
and set C = 0
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (right_inlet_bdy, uv, &inflow_right));
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (right_inlet_bdy, C_only, &zero));
Note that the remaining boundary conditions are the natural boundary conditions for the concentration
on the wall (grad(c) dot n = 0) and natural boundary conditions for the velocity and the concentration
on the outlets ((grad(velocity) dot n - p n) dot t = 0, grad(C) dot n = 0)
Do the parent's initialization after variables and boundary constraints are defined
Do the parent's initialization after variables and boundary constraints are defined
FEMSystem::init_data();
}
void CoupledSystem::init_context(DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
We should prerequest all the data
we will need to build the linear system.
Note that the concentration and velocity components
use the same basis.
c.element_fe_var[u_var]->get_JxW();
c.element_fe_var[u_var]->get_phi();
c.element_fe_var[u_var]->get_dphi();
c.element_fe_var[u_var]->get_xyz();
c.element_fe_var[p_var]->get_phi();
c.side_fe_var[u_var]->get_JxW();
c.side_fe_var[u_var]->get_phi();
c.side_fe_var[u_var]->get_xyz();
}
bool CoupledSystem::element_time_derivative (bool request_jacobian,
DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
First we get some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weights for interior integration
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW =
c.element_fe_var[u_var]->get_JxW();
The velocity shape functions at interior quadrature points.
const std::vector<std::vector<Real> >& phi =
c.element_fe_var[u_var]->get_phi();
The velocity shape function gradients at interior
quadrature points.
const std::vector<std::vector<RealGradient> >& dphi =
c.element_fe_var[u_var]->get_dphi();
The pressure shape functions at interior
quadrature points.
const std::vector<std::vector<Real> >& psi =
c.element_fe_var[p_var]->get_phi();
The number of local degrees of freedom in each variable
const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
libmesh_assert_equal_to (n_u_dofs, c.dof_indices_var[v_var].size());
The subvectors and submatrices we need to fill:
DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
DenseSubMatrix<Number> &Kup = *c.elem_subjacobians[u_var][p_var];
DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
DenseSubMatrix<Number> &Kvv = *c.elem_subjacobians[v_var][v_var];
DenseSubMatrix<Number> &Kvp = *c.elem_subjacobians[v_var][p_var];
DenseSubVector<Number> &Fv = *c.elem_subresiduals[v_var];
DenseSubMatrix<Number> &KCu = *c.elem_subjacobians[C_var][u_var];
DenseSubMatrix<Number> &KCv = *c.elem_subjacobians[C_var][v_var];
DenseSubMatrix<Number> &KCC = *c.elem_subjacobians[C_var][C_var];
DenseSubVector<Number> &FC = *c.elem_subresiduals[C_var];
Now we will build the element Jacobian and residual.
Constructing the residual requires the solution and its
gradient from the previous timestep. This must be
calculated at each quadrature point by summing the
solution degree-of-freedom values by the appropriate
weight functions.
unsigned int n_qpoints = c.element_qrule->n_points();
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Compute the solution & its gradient at the old Newton iterate
Number p = c.interior_value(p_var, qp),
u = c.interior_value(u_var, qp),
v = c.interior_value(v_var, qp);
Gradient grad_u = c.interior_gradient(u_var, qp),
grad_v = c.interior_gradient(v_var, qp),
grad_C = c.interior_gradient(C_var, qp);
Definitions for convenience. It is sometimes simpler to do a
dot product if you have the full vector at your disposal.
NumberVectorValue U (u, v);
const Number C_x = grad_C(0);
const Number C_y = grad_C(1);
First, an i-loop over the velocity degrees of freedom.
We know that n_u_dofs == n_v_dofs so we can compute contributions
for both at the same time.
for (unsigned int i=0; i != n_u_dofs; i++)
{
Stokes equations residuals
Fu(i) += JxW[qp] *
(p*dphi[i][qp](0) - // pressure term
(grad_u*dphi[i][qp])); // diffusion term
Fv(i) += JxW[qp] *
(p*dphi[i][qp](1) - // pressure term
(grad_v*dphi[i][qp])); // diffusion term
Concentration Equation Residual
FC(i) += JxW[qp] *
( (U*grad_C)*phi[i][qp] + // convection term
(1./Peclet)*(grad_C*dphi[i][qp]) ); // diffusion term
Note that the Fp block is identically zero unless we are using
some kind of artificial compressibility scheme...
if (request_jacobian && c.elem_solution_derivative)
{
libmesh_assert (c.elem_solution_derivative == 1.0);
Matrix contributions for the uu and vv couplings.
for (unsigned int j=0; j != n_u_dofs; j++)
{
Kuu(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term */
Kvv(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term */
KCu(i,j) += JxW[qp]* ( (phi[j][qp]*C_x)*phi[i][qp] ); /* convection term */
KCv(i,j) += JxW[qp]*( (phi[j][qp]*C_y)*phi[i][qp] ); /* convection term */
KCC(i,j) += JxW[qp]*
( (U*dphi[j][qp])*phi[i][qp] + /* nonlinear term (convection) */
(1./Peclet)*(dphi[j][qp]*dphi[i][qp]) ); /* diffusion term */
}
Matrix contributions for the up and vp couplings.
for (unsigned int j=0; j != n_p_dofs; j++)
{
Kup(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](0);
Kvp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](1);
}
}
}
} // end of the quadrature point qp-loop
return request_jacobian;
}
bool CoupledSystem::element_constraint (bool request_jacobian,
DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
Here we define some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weight for interior integration
Element Jacobian * quadrature weight for interior integration
const std::vector<Real> &JxW = c.element_fe_var[u_var]->get_JxW();
The velocity shape function gradients at interior
quadrature points.
const std::vector<std::vector<RealGradient> >& dphi =
c.element_fe_var[u_var]->get_dphi();
The pressure shape functions at interior
quadrature points.
const std::vector<std::vector<Real> >& psi =
c.element_fe_var[p_var]->get_phi();
The number of local degrees of freedom in each variable
const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
The subvectors and submatrices we need to fill:
DenseSubMatrix<Number> &Kpu = *c.elem_subjacobians[p_var][u_var];
DenseSubMatrix<Number> &Kpv = *c.elem_subjacobians[p_var][v_var];
DenseSubVector<Number> &Fp = *c.elem_subresiduals[p_var];
Add the constraint given by the continuity equation
unsigned int n_qpoints = c.element_qrule->n_points();
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Compute the velocity gradient at the old Newton iterate
Gradient grad_u = c.interior_gradient(u_var, qp),
grad_v = c.interior_gradient(v_var, qp);
Now a loop over the pressure degrees of freedom. This
computes the contributions of the continuity equation.
for (unsigned int i=0; i != n_p_dofs; i++)
{
Fp(i) += JxW[qp] * psi[i][qp] *
(grad_u(0) + grad_v(1));
if (request_jacobian && c.elem_solution_derivative)
{
libmesh_assert (c.elem_solution_derivative == 1.0);
for (unsigned int j=0; j != n_u_dofs; j++)
{
Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
}
}
}
} // end of the quadrature point qp-loop
return request_jacobian;
}
void CoupledSystem::postprocess()
{
We need to overload the postprocess function to set the computed_QoI variable of the CoupledSystem class
to the qoi value stored in System::qoi[0]
computed_QoI = 0.0;
computed_QoI = System::qoi[0];
}
These functions supply the nonlinear weighting for the adjoint residual error estimate which
arise due to the convection term in the convection-diffusion equation:
||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} + ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2}
||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} + ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2}
These functions compute (u_1)_h or C,1_h , and (u_2)_h or C,2_h , and supply it to the weighted patch recovery error estimator
In CoupledFEMFunctionsx, the object built with var = 0, returns the (u_1)_h weight, while
the object built with var = 1, returns the C,1_h weight. The switch statment
distinguishes the behavior of the two objects
Same thing for CoupledFEMFunctionsy
Number CoupledFEMFunctionsx::operator()(const FEMContext& c, const Point& p,
const Real /* time */)
{
Number weight = 0.0;
switch(var)
{
case 0:
{
Gradient grad_C = c.point_gradient(3, p);
weight = grad_C(0);
}
break;
case 3:
{
Number u = c.point_value(0, p);
weight = u;
}
break;
default:
{
std::cout<<"Wrong variable number"<<var<<" passed to CoupledFEMFunctionsx object ! Quitting !"<<std::endl;
libmesh_error();
}
}
return weight;
}
Number CoupledFEMFunctionsy::operator()(const FEMContext& c, const Point& p,
const Real /* time */)
{
Number weight = 0.0;
switch(var)
{
case 1:
{
Gradient grad_C = c.point_gradient(3, p);
weight = grad_C(1);
}
break;
case 3:
{
Number v = c.point_value(1, p);
weight = v;
}
break;
default:
{
std::cout<<"Wrong variable number "<<var<<" passed to CoupledFEMFunctionsy object ! Quitting !"<<std::endl;
libmesh_error();
}
}
return weight;
}
The source file domain.C with comments:
#include "libmesh/mesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_refinement.h"
#include "domain.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Mesh construction
void build_domain (Mesh &mesh, FEMParameters ¶m)
{
mesh.read(param.domainfile);
std::cout<<"Making elements 2nd order"<<std::endl;
Right now we are setting approximation orders in the code, rather than reading them in
That needs to be fixed and the second ordering should be done only if one of the
approximation orders is greater than 1
mesh.all_second_order();
}
The source file femparameters.C with comments:
#include "femparameters.h"
using namespace libMesh;
#define GETPOT_INPUT(A) { A = input(#A, A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_INT_INPUT(A) { A = input(#A, (int)A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_REGISTER(A) { \
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
void FEMParameters::read(GetPot &input)
{
std::vector<std::string> variable_names;
GETPOT_INT_INPUT(initial_timestep);
GETPOT_INT_INPUT(n_timesteps);
GETPOT_INPUT(transient);
GETPOT_INT_INPUT(deltat_reductions);
GETPOT_INPUT(timesolver_core);
GETPOT_INPUT(end_time);
GETPOT_INPUT(deltat);
GETPOT_INPUT(timesolver_theta);
GETPOT_INPUT(timesolver_maxgrowth);
GETPOT_INPUT(timesolver_upper_tolerance);
GETPOT_INPUT(timesolver_tolerance);
GETPOT_INPUT(steadystate_tolerance);
GETPOT_REGISTER(timesolver_norm);
const unsigned int n_timesolver_norm = input.vector_variable_size("timesolver_norm");
timesolver_norm.resize(n_timesolver_norm, L2);
for (unsigned int i=0; i != n_timesolver_norm; ++i)
{
int current_norm = 0; // L2
if (timesolver_norm[i] == H1)
current_norm = 1;
if (timesolver_norm[i] == H2)
current_norm = 2;
current_norm = input("timesolver_norm", current_norm, i);
if (current_norm == 0)
timesolver_norm[i] = L2;
else if (current_norm == 1)
timesolver_norm[i] = H1;
else if (current_norm == 2)
timesolver_norm[i] = H2;
else
timesolver_norm[i] = DISCRETE_L2;
}
GETPOT_INT_INPUT(dimension);
GETPOT_INPUT(domaintype);
GETPOT_INPUT(domainfile);
GETPOT_INPUT(elementtype);
GETPOT_INPUT(fine_mesh_file_primal);
GETPOT_INPUT(fine_mesh_soln_primal);
GETPOT_INPUT(fine_mesh_file_adjoint);
GETPOT_INPUT(fine_mesh_soln_adjoint);
GETPOT_INPUT(elementorder);
GETPOT_INPUT(domain_xmin);
GETPOT_INPUT(domain_ymin);
GETPOT_INPUT(domain_zmin);
GETPOT_INPUT(domain_edge_width);
GETPOT_INPUT(domain_edge_length);
GETPOT_INPUT(domain_edge_height);
GETPOT_INT_INPUT(coarsegridx);
GETPOT_INT_INPUT(coarsegridy);
GETPOT_INT_INPUT(coarsegridz);
GETPOT_INT_INPUT(coarserefinements);
GETPOT_INT_INPUT(coarsecoarsenings);
GETPOT_INT_INPUT(extrarefinements);
GETPOT_INPUT(use_petsc_snes);
GETPOT_INPUT(time_solver_quiet);
GETPOT_INPUT(solver_quiet);
GETPOT_INPUT(reuse_preconditioner);
GETPOT_INPUT(require_residual_reduction);
GETPOT_INPUT(min_step_length);
GETPOT_INT_INPUT(max_linear_iterations);
GETPOT_INT_INPUT(max_nonlinear_iterations);
GETPOT_INPUT(relative_step_tolerance);
GETPOT_INPUT(relative_residual_tolerance);
GETPOT_INPUT(initial_linear_tolerance);
GETPOT_INPUT(minimum_linear_tolerance);
GETPOT_INPUT(linear_tolerance_multiplier);
GETPOT_INT_INPUT(nelem_target);
GETPOT_INPUT(global_tolerance);
GETPOT_INPUT(refine_fraction);
GETPOT_INPUT(coarsen_fraction);
GETPOT_INPUT(coarsen_threshold);
GETPOT_INT_INPUT(max_adaptivesteps);
GETPOT_INT_INPUT(initial_adaptivesteps);
GETPOT_INT_INPUT(initial_sobolev_order);
GETPOT_INT_INPUT(initial_extra_quadrature);
GETPOT_INPUT(refine_uniformly);
GETPOT_INPUT(indicator_type);
GETPOT_INPUT(adjoint_residual_type);
GETPOT_INPUT(patch_reuse);
GETPOT_INT_INPUT(sobolev_order);
GETPOT_INPUT(alternate_with_uniform_steps);
GETPOT_INPUT(alternate_step_number);
GETPOT_INPUT(component_wise_error);
GETPOT_INPUT(compute_sensitivities);
GETPOT_INPUT(compare_to_fine_solution_primal);
GETPOT_INPUT(compare_to_fine_solution_adjoint);
GETPOT_INPUT(do_forward_sensitivity);
GETPOT_INPUT(do_adjoint_sensitivity);
GETPOT_INPUT(postprocess_adjoint);
GETPOT_INT_INPUT(write_interval);
GETPOT_INPUT(output_xda);
GETPOT_INPUT(output_xdr);
GETPOT_INPUT(output_gz);
#ifndef LIBMESH_HAVE_GZSTREAM
output_gz = false;
#endif
GETPOT_INPUT(output_bz2);
#ifndef LIBMESH_HAVE_BZ2
output_bz2 = false;
#endif
GETPOT_INPUT(output_gmv);
GETPOT_INPUT(write_gmv_error);
#ifndef LIBMESH_HAVE_GMV
output_gmv = false;
write_gmv_error = false;
#endif
GETPOT_INPUT(output_tecplot);
GETPOT_INPUT(write_tecplot_error);
#ifndef LIBMESH_HAVE_TECPLOT_API
output_tecplot = false;
write_tecplot_error = false;
#endif
GETPOT_INPUT(run_simulation);
GETPOT_INPUT(run_postprocess);
run_simulation = input("run_simulation", run_simulation);
run_postprocess = input("run_postprocess", run_postprocess);
GETPOT_REGISTER(fe_family);
const unsigned int n_fe_family =
std::max(1u, input.vector_variable_size("fe_family"));
fe_family.resize(n_fe_family, "LAGRANGE");
for (unsigned int i=0; i != n_fe_family; ++i)
fe_family[i] = input("fe_family", fe_family[i].c_str(), i);
GETPOT_REGISTER(fe_order);
const unsigned int n_fe_order =
input.vector_variable_size("fe_order");
fe_order.resize(n_fe_order, 1);
for (unsigned int i=0; i != n_fe_order; ++i)
fe_order[i] = input("fe_order", (int)fe_order[i], i);
GETPOT_INPUT(extra_quadrature_order);
GETPOT_INPUT(analytic_jacobians);
GETPOT_INPUT(verify_analytic_jacobians);
GETPOT_INPUT(print_solution_norms);
GETPOT_INPUT(print_solutions);
GETPOT_INPUT(print_residual_norms);
GETPOT_INPUT(print_residuals);
GETPOT_INPUT(print_jacobian_norms);
GETPOT_INPUT(print_jacobians);
std::vector bad_variables =
input.unidentified_arguments(variable_names);
if (libMesh::processor_id() == 0 && !bad_variables.empty()) { std::cerr << "ERROR: Unrecognized variables:" << std::endl; for (unsigned int i = 0; i != bad_variables.size(); ++i) std::cerr << bad_variables[i] << std::endl; std::cerr << "not found among recognized variables." << std::endl; for (unsigned int i = 0; i != variable_names.size(); ++i) std::cerr << variable_names[i] << std::endl; libmesh_error(); }
if (libMesh::processor_id() == 0 && !bad_variables.empty()) { std::cerr << "ERROR: Unrecognized variables:" << std::endl; for (unsigned int i = 0; i != bad_variables.size(); ++i) std::cerr << bad_variables[i] << std::endl; std::cerr << "not found among recognized variables." << std::endl; for (unsigned int i = 0; i != variable_names.size(); ++i) std::cerr << variable_names[i] << std::endl; libmesh_error(); }
}
The source file H-qoi.C with comments:
#include "H-qoi.h"
Here we define the functions to compute the QoI (side_qoi)
and supply the right hand side for the associated adjoint problem (side_qoi_derivative)
using namespace libMesh;
void CoupledSystemQoI::init_qoi( std::vector<Number>& sys_qoi)
{
Only 1 qoi to worry about
sys_qoi.resize(1);
return;
}
This function supplies the right hand side for the adjoint problem
We only have one QoI, so we don't bother checking the qois argument
to see if it was requested from us
void CoupledSystemQoI::side_qoi_derivative (DiffContext &context,
const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
Get velocity basis functions phi
const std::vector<std::vector<Real> > &phi = c.side_fe_var[0]->get_phi();
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
The number of local degrees of freedom in each variable
const unsigned int n_u_dofs = c.dof_indices_var[1].size();
DenseSubVector<Number> &Qu = *c.elem_qoi_subderivatives[0][0];
DenseSubVector<Number> &QC = *c.elem_qoi_subderivatives[0][3];
Now we will build the element Jacobian and residual.
Constructing the residual requires the solution and its
gradient from the previous timestep. This must be
calculated at each quadrature point by summing the
solution degree-of-freedom values by the appropriate
weight functions.
unsigned int n_qpoints = c.side_qrule->n_points();
Number u = 0. ;
Number C = 0. ;
// If side is on outlet
if(c.has_side_boundary_id(2))
{
Loop over all the qps on this side
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Real x = q_point[qp](0);
If side is on left outlet
if(x < 0.)
{
Get u at the qp
u = c.side_value(0,qp);
C = c.side_value(3,qp);
Add the contribution from each basis function
for (unsigned int i=0; i != n_u_dofs; i++)
{
Qu(i) += JxW[qp] * -phi[i][qp] * C;
QC(i) += JxW[qp] * phi[i][qp] * -u;
}
} // end if
} // end quadrature loop
} // end if on outlet
}
This function computes the actual QoI
void CoupledSystemQoI::side_qoi(DiffContext &context, const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
First we get some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weights for interior integration
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
Loop over qp's, compute the function at each qp and add
to get the QoI
unsigned int n_qpoints = c.side_qrule->n_points();
Number dQoI_0 = 0. ;
Number u = 0. ;
Number C = 0. ;
If side is on the left outlet
if(c.has_side_boundary_id(2))
{
Loop over all the qps on this side
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Real x = q_point[qp](0);
if(x < 0.)
{
Get u and C at the qp
u = c.side_value(0,qp);
C = c.side_value(3,qp);
dQoI_0 += JxW[qp] * -u * C;
} // end if
} // end quadrature loop
} // end if on bdry
c.elem_qoi[0] += dQoI_0;
}
The source file initial.C with comments:
#include "initial.h"
#include "femparameters.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
void read_initial_parameters()
{
}
void finish_initialization()
{
}
Initial conditions
Number initial_value(const Point& /* p */,
const Parameters& /* param */,
const std::string&,
const std::string&)
{
return 1.;
}
Gradient initial_grad(const Point& /* p */,
const Parameters& /* param */,
const std::string&,
const std::string&)
{
return 0.;
}
The source file output.C with comments:
#include <fstream>
#include <unistd.h>
#include "libmesh/libmesh_common.h"
#include "output.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
void start_output(unsigned int timesteps,
std::string filename,
std::string varname)
{
Truncate then append if we're doing a restart
if (timesteps)
{
std::ifstream infile (filename.c_str());
char buf[1025];
if (!infile.is_open())
libmesh_error();
Count off the lines we want to save, including
the original header
for (unsigned int i=0; i != timesteps+1; ++i)
infile.getline(buf, 1024);
if (!infile.good())
libmesh_error();
unsigned int length = infile.tellg();
infile.close();
Then throw away the rest
int err = truncate(filename.c_str(), length);
if (err != 0)
libmesh_error();
}
Write new headers if we're starting a new simulation
else
{
std::ofstream outfile (filename.c_str(), std::ios_base::out);
outfile << varname << " = [" << std::endl;
}
}
The source file coupled_system.h without comments:
#ifndef __coupled_system_h__
#define __coupled_system_h__
#include "libmesh/fem_function_base.h"
#include "libmesh/fem_system.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/parameter_vector.h"
using namespace libMesh;
class CoupledSystem : public FEMSystem
{
public:
CoupledSystem(EquationSystems& es,
const std::string& name_in,
const unsigned int number_in)
: FEMSystem(es, name_in, number_in), Peclet(1.) {qoi.resize(1);}
Number &get_QoI_value()
{
return computed_QoI;
}
Number &get_parameter_value(unsigned int parameter_index)
{
return parameters[parameter_index];
}
ParameterVector &get_parameter_vector()
{
parameter_vector.resize(parameters.size());
for(unsigned int i = 0; i != parameters.size(); ++i)
{
parameter_vector[i] = ¶meters[i];
}
return parameter_vector;
}
Real &get_Pe()
{
return Peclet;
}
protected:
virtual void init_data ();
virtual void init_context(DiffContext &context);
virtual bool element_time_derivative (bool request_jacobian,
DiffContext& context);
virtual bool element_constraint (bool request_jacobian,
DiffContext& context);
virtual void postprocess ();
std::vector<Number> parameters;
unsigned int p_var, u_var, v_var, C_var;
ParameterVector parameter_vector;
Real Peclet;
Number computed_QoI;
};
class CoupledFEMFunctionsx : public FEMFunctionBase<Number>
{
public:
CoupledFEMFunctionsx(System& /* sys */, unsigned int var_number) {var = var_number;}
virtual ~CoupledFEMFunctionsx () {}
virtual AutoPtr<FEMFunctionBase<Number> > clone () const
{return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsx(*this) ); }
virtual void operator() (const FEMContext&, const Point&,
const Real, DenseVector<Number>&)
{libmesh_error();}
virtual Number operator() (const FEMContext&, const Point& p,
const Real time = 0.);
private:
unsigned int var;
};
class CoupledFEMFunctionsy : public FEMFunctionBase<Number>
{
public:
CoupledFEMFunctionsy(System& /* sys */, unsigned int var_number) {var = var_number;}
virtual ~CoupledFEMFunctionsy () {}
virtual AutoPtr<FEMFunctionBase<Number> > clone () const
{return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsy(*this) ); }
virtual void operator() (const FEMContext&, const Point&,
const Real,
DenseVector<Number>&)
{libmesh_error();}
virtual Number operator() (const FEMContext&, const Point& p,
const Real time = 0.);
private:
unsigned int var;
};
#endif //__coupled_system_h__
The source file domain.h without comments:
#include "femparameters.h" #include "libmesh/mesh.h" class FEMParameters; void build_domain (Mesh &mesh, FEMParameters ¶m);
The source file femparameters.h without comments:
#ifndef __fem_parameters_h__
#define __fem_parameters_h__
#include <limits>
#include "libmesh/libmesh_common.h"
#include "libmesh/dof_map.h"
#include "libmesh/enum_norm_type.h"
#include "libmesh/getpot.h"
using namespace libMesh;
class FEMParameters
{
public:
FEMParameters() :
initial_timestep(0), n_timesteps(100),
transient(true),
deltat_reductions(0),
timesolver_core("euler"),
end_time(std::numeric_limits<Real>::max()),
deltat(0.0001), timesolver_theta(0.5),
timesolver_maxgrowth(0.), timesolver_tolerance(0.),
timesolver_upper_tolerance(0.),
steadystate_tolerance(0.),
timesolver_norm(0,L2),
dimension(2),
domaintype("square"), domainfile("mesh.xda"), elementtype("quad"),
fine_mesh_file_primal("fine_mesh.xda"), fine_mesh_soln_primal("fine_mesh_soln.xda"),
fine_mesh_file_adjoint("fine_mesh.xda"), fine_mesh_soln_adjoint("fine_mesh_soln.xda"),
elementorder(2),
domain_xmin(0.0), domain_ymin(0.0), domain_zmin(0.0),
domain_edge_width(1.0), domain_edge_length(1.0), domain_edge_height(1.0),
coarsegridx(1), coarsegridy(1), coarsegridz(1),
coarserefinements(0), coarsecoarsenings(0), extrarefinements(0),
use_petsc_snes(false),
time_solver_quiet(true), solver_quiet(true),
reuse_preconditioner(false),
require_residual_reduction(true),
min_step_length(1e-5),
max_linear_iterations(200000), max_nonlinear_iterations(20),
relative_step_tolerance(1.e-7), relative_residual_tolerance(1.e-10),
initial_linear_tolerance(1.e-3), minimum_linear_tolerance(TOLERANCE*TOLERANCE),
linear_tolerance_multiplier(1.e-3),
nelem_target(30000), global_tolerance(0.0),
refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
refine_uniformly(false),
max_adaptivesteps(1),
initial_adaptivesteps(0), initial_sobolev_order(1),
initial_extra_quadrature(0),
indicator_type("kelly"), patch_reuse(true), sobolev_order(1), adjoint_residual_type("patchpatch"), alternate_with_uniform_steps("false"), alternate_step_number(2),
component_wise_error("false"), compare_to_fine_solution_primal(false), compare_to_fine_solution_adjoint(false), compute_sensitivities(false), do_forward_sensitivity(true), do_adjoint_sensitivity(false),
postprocess_adjoint(false),
write_interval(10),
write_gmv_error(false), write_tecplot_error(false),
output_xda(false), output_xdr(false),
output_bz2(true), output_gz(true),
output_gmv(false), output_tecplot(false),
run_simulation(true), run_postprocess(false),
fe_family(1, "LAGRANGE"), fe_order(1, 1),
extra_quadrature_order(0),
analytic_jacobians(true), verify_analytic_jacobians(0.0),
print_solution_norms(false), print_solutions(false),
print_residual_norms(false), print_residuals(false),
print_jacobian_norms(false), print_jacobians(false) {}
void read(GetPot &input);
unsigned int initial_timestep, n_timesteps;
bool transient;
unsigned int deltat_reductions;
std::string timesolver_core;
Real end_time, deltat, timesolver_theta,
timesolver_maxgrowth, timesolver_tolerance,
timesolver_upper_tolerance, steadystate_tolerance;
std::vector<FEMNormType> timesolver_norm;
unsigned int dimension;
std::string domaintype, domainfile, elementtype, fine_mesh_file_primal, fine_mesh_soln_primal, fine_mesh_file_adjoint, fine_mesh_soln_adjoint;
Real elementorder;
Real domain_xmin, domain_ymin, domain_zmin;
Real domain_edge_width, domain_edge_length, domain_edge_height;
unsigned int coarsegridx, coarsegridy, coarsegridz;
unsigned int coarserefinements, coarsecoarsenings, extrarefinements;
bool use_petsc_snes;
bool time_solver_quiet, solver_quiet, reuse_preconditioner, require_residual_reduction;
Real min_step_length;
unsigned int max_linear_iterations, max_nonlinear_iterations;
Real relative_step_tolerance, relative_residual_tolerance,
initial_linear_tolerance, minimum_linear_tolerance,
linear_tolerance_multiplier;
unsigned int nelem_target;
Real global_tolerance;
Real refine_fraction, coarsen_fraction, coarsen_threshold;
bool refine_uniformly;
unsigned int max_adaptivesteps;
unsigned int initial_adaptivesteps;
unsigned int initial_sobolev_order;
unsigned int initial_extra_quadrature;
std::string indicator_type;
bool patch_reuse;
unsigned int sobolev_order;
std::string adjoint_residual_type;
bool alternate_with_uniform_steps;
unsigned int alternate_step_number;
bool component_wise_error;
bool compare_to_fine_solution_primal, compare_to_fine_solution_adjoint;
bool compute_sensitivities;
bool do_forward_sensitivity;
bool do_adjoint_sensitivity;
bool postprocess_adjoint;
unsigned int write_interval;
bool write_gmv_error, write_tecplot_error, output_xda, output_xdr,
output_bz2, output_gz, output_gmv, output_tecplot;
bool run_simulation, run_postprocess;
std::vector<std::string> fe_family;
std::vector<unsigned int> fe_order;
int extra_quadrature_order;
bool analytic_jacobians;
Real verify_analytic_jacobians;
bool print_solution_norms, print_solutions,
print_residual_norms, print_residuals,
print_jacobian_norms, print_jacobians;
};
#endif // __fem_parameters_h__
The source file H-qoi.h without comments:
#ifndef H_QOI_H
#define H_QOI_H
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
#include "libmesh/diff_qoi.h"
using namespace libMesh;
class CoupledSystemQoI : public DifferentiableQoI
{
public:
CoupledSystemQoI(){}
virtual ~CoupledSystemQoI(){}
virtual void init_qoi( std::vector<Number>& sys_qoi);
virtual void postprocess( ){}
virtual void side_qoi_derivative(DiffContext &context, const QoISet & qois);
virtual void side_qoi(DiffContext &context, const QoISet & qois);
virtual AutoPtr<DifferentiableQoI> clone( )
{ AutoPtr<DifferentiableQoI> my_clone( new CoupledSystemQoI );
*my_clone = *this;
return my_clone;
}
};
#endif // H_QOI_H
The source file initial.h without comments:
#include "libmesh/parameters.h"
#include "libmesh/point.h"
#include "libmesh/vector_value.h"
using namespace libMesh;
void read_initial_parameters();
void finish_initialization();
Number initial_value(const Point& p,
const Parameters&,
const std::string&,
const std::string&);
Gradient initial_grad(const Point& p,
const Parameters&,
const std::string&,
const std::string&);
The source file mysystems.h without comments:
#include "coupled_system.h"
#include "femparameters.h"
using namespace libMesh;
FEMSystem &build_system(EquationSystems &es, GetPot &, FEMParameters& /* param */)
{
CoupledSystem &system = es.add_system<CoupledSystem> ("CoupledSystem");
return system;
}
Number exact_value(const Point& /* p */, // xyz location
const Parameters& /* param */, // EquationSystem parameters
const std::string&, // sys_name
const std::string&) // unknown_name
{
std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
return 0;
}
Gradient exact_grad(const Point& /* p */, // xyz location
const Parameters& /* param */, // EquationSystems parameters
const std::string&, // sys_name
const std::string&) // unknown_name
{
std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
return 0;
}
The source file output.h without comments:
#include <string>
using namespace libMesh;
void start_output(unsigned int timesteps,
std::string filename,
std::string varname);
The source file adjoints_ex3.C without comments:
#include <iostream>
#include <sys/time.h>
#include <iomanip>
#include "libmesh/equation_systems.h"
#include "libmesh/twostep_time_solver.h"
#include "libmesh/euler_solver.h"
#include "libmesh/euler2_solver.h"
#include "libmesh/steady_solver.h"
#include "libmesh/newton_solver.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/petsc_diff_solver.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_tools.h"
#include "libmesh/mesh_base.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/system.h"
#include "libmesh/system_norm.h"
#include "libmesh/adjoint_residual_error_estimator.h"
#include "libmesh/const_fem_function.h"
#include "libmesh/error_vector.h"
#include "libmesh/fem_function_base.h"
#include "libmesh/getpot.h"
#include "libmesh/gmv_io.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/parameter_vector.h"
#include "libmesh/patch_recovery_error_estimator.h"
#include "libmesh/petsc_vector.h"
#include "libmesh/sensitivity_data.h"
#include "libmesh/tecplot_io.h"
#include "libmesh/uniform_refinement_estimator.h"
#include "libmesh/qoi_set.h"
#include "libmesh/weighted_patch_recovery_error_estimator.h"
#include "coupled_system.h"
#include "domain.h"
#include "initial.h"
#include "femparameters.h"
#include "mysystems.h"
#include "output.h"
#include "H-qoi.h"
using namespace libMesh;
std::string numbered_filename(unsigned int t_step, // The timestep count
unsigned int a_step, // The adaptive step count
std::string solution_type, // primal or adjoint solve
std::string type,
std::string extension,
FEMParameters ¶m)
{
std::ostringstream file_name;
file_name << solution_type
<< ".out."
<< type
<< '.'
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step
<< '.'
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step
<< '.'
<< extension;
if (param.output_bz2)
file_name << ".bz2";
else if (param.output_gz)
file_name << ".gz";
return file_name.str();
}
void write_output(EquationSystems &es,
unsigned int t_step, // The timestep count
unsigned int a_step, // The adaptive step count
std::string solution_type, // primal or adjoint solve
FEMParameters ¶m)
{
MeshBase &mesh = es.get_mesh();
#ifdef LIBMESH_HAVE_GMV
if (param.output_gmv)
{
std::ostringstream file_name_gmv;
file_name_gmv << solution_type
<< ".out.gmv."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step
<< '.'
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step;
GMVIO(mesh).write_equation_systems
(file_name_gmv.str(), es);
}
#endif
#ifdef LIBMESH_HAVE_TECPLOT_API
if (param.output_tecplot)
{
std::ostringstream file_name_tecplot;
file_name_tecplot << solution_type
<< ".out."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step
<< '.'
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step
<< ".plt";
TecplotIO(mesh).write_equation_systems
(file_name_tecplot.str(), es);
}
#endif
if (param.output_xda || param.output_xdr)
mesh.renumber_nodes_and_elements();
if (param.output_xda)
{
mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param));
es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xda", param),
libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
EquationSystems::WRITE_ADDITIONAL_DATA);
}
if (param.output_xdr)
{
mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param));
es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param),
libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
EquationSystems::WRITE_ADDITIONAL_DATA);
}
}
void write_output_headers(FEMParameters ¶m)
{
if (libMesh::processor_id() != 0)
return;
start_output(param.initial_timestep, "out_clocktime.m", "vector_clocktime");
if (param.run_simulation)
{
start_output(param.initial_timestep, "out_activemesh.m", "table_activemesh");
start_output(param.initial_timestep, "out_solvesteps.m", "table_solvesteps");
if (param.timesolver_tolerance)
{
start_output(param.initial_timestep, "out_time.m", "vector_time");
start_output(param.initial_timestep, "out_timesteps.m", "vector_timesteps");
}
if (param.steadystate_tolerance)
start_output(param.initial_timestep, "out_changerate.m", "vector_changerate");
}
}
void write_output_solvedata(EquationSystems &es,
unsigned int a_step,
unsigned int newton_steps,
unsigned int krylov_steps,
unsigned int tv_sec,
unsigned int tv_usec)
{
MeshBase &mesh = es.get_mesh();
unsigned int n_active_elem = mesh.n_active_elem();
unsigned int n_active_dofs = es.n_active_dofs();
if (libMesh::processor_id() == 0)
{
std::ofstream activemesh ("out_activemesh.m",
std::ios_base::app | std::ios_base::out);
activemesh.precision(17);
activemesh << (a_step + 1) << ' '
<< n_active_elem << ' '
<< n_active_dofs << std::endl;
std::ofstream solvesteps ("out_solvesteps.m",
std::ios_base::app | std::ios_base::out);
solvesteps.precision(17);
solvesteps << newton_steps << ' '
<< krylov_steps << std::endl;
std::ofstream clocktime ("out_clocktime.m",
std::ios_base::app | std::ios_base::out);
clocktime.precision(17);
clocktime << tv_sec << '.' << tv_usec << std::endl;
}
}
void write_output_footers(FEMParameters ¶m)
{
if (libMesh::processor_id() == 0)
{
std::ofstream clocktime ("out_clocktime.m",
std::ios_base::app | std::ios_base::out);
clocktime << "];" << std::endl;
if (param.run_simulation)
{
std::ofstream activemesh ("out_activemesh.m",
std::ios_base::app | std::ios_base::out);
activemesh << "];" << std::endl;
std::ofstream solvesteps ("out_solvesteps.m",
std::ios_base::app | std::ios_base::out);
solvesteps << "];" << std::endl;
if (param.timesolver_tolerance)
{
std::ofstream times ("out_time.m",
std::ios_base::app | std::ios_base::out);
times << "];" << std::endl;
std::ofstream timesteps ("out_timesteps.m",
std::ios_base::app | std::ios_base::out);
timesteps << "];" << std::endl;
}
if (param.steadystate_tolerance)
{
std::ofstream changerate ("out_changerate.m",
std::ios_base::app | std::ios_base::out);
changerate << "];" << std::endl;
}
}
std::ofstream complete ("complete");
complete << "complete" << std::endl;
}
}
#if defined(LIBMESH_HAVE_GMV) || defined(LIBMESH_HAVE_TECPLOT_API)
void write_error(EquationSystems &es,
ErrorVector &error,
unsigned int t_number,
unsigned int a_number,
FEMParameters ¶m,
std::string error_type)
#else
void write_error(EquationSystems &,
ErrorVector &,
unsigned int,
unsigned int,
FEMParameters &,
std::string)
#endif
{
#ifdef LIBMESH_HAVE_GMV
if (param.write_gmv_error)
{
std::ostringstream error_gmv;
error_gmv << "error.gmv."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< a_number
<< "."
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< t_number
<< error_type;
error.plot_error(error_gmv.str(), es.get_mesh());
}
#endif
#ifdef LIBMESH_HAVE_TECPLOT_API
if (param.write_tecplot_error)
{
std::ostringstream error_tecplot;
error_tecplot << "error.plt."
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< a_number
<< "."
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< t_number
<< error_type;
error.plot_error(error_tecplot.str(), es.get_mesh());
}
#endif
}
void read_output(EquationSystems &es,
unsigned int t_step,
unsigned int a_step,
std::string solution_type,
FEMParameters ¶m)
{
MeshBase &mesh = es.get_mesh();
std::string file_name_mesh, file_name_soln;
if (param.output_xda)
{
file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param);
file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xda", param);
}
else if (param.output_xdr)
{
file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param);
file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param);
}
mesh.read(file_name_mesh);
es.read(file_name_soln, libMeshEnums::READ,
EquationSystems::READ_HEADER |
EquationSystems::READ_DATA |
EquationSystems::READ_ADDITIONAL_DATA);
for (unsigned int i = 0; i != es.n_systems(); ++i)
es.get_system<FEMSystem>(i).update();
Real current_time = 0., current_timestep = 0.;
if (param.timesolver_tolerance)
{
std::ifstream times ("out_time.m");
std::ifstream timesteps ("out_timesteps.m");
if (times.is_open() && timesteps.is_open())
{
const unsigned int headersize = 25;
char header[headersize];
timesteps.getline (header, headersize);
if (strcmp(header, "vector_timesteps = [") != 0)
{
std::cout << "Bad header in out_timesteps.m:" << std::endl
<< header
<< std::endl;
libmesh_error();
}
times.getline (header, headersize);
if (strcmp(header, "vector_time = [") != 0)
{
std::cout << "Bad header in out_time.m:" << std::endl
<< header
<< std::endl;
libmesh_error();
}
for (unsigned int i = 0; i != t_step; ++i)
{
if (!times.good())
libmesh_error();
times >> current_time;
timesteps >> current_timestep;
}
current_time += current_timestep;
}
else
libmesh_error();
}
else
current_time = t_step * param.deltat;
for (unsigned int i = 0; i != es.n_systems(); ++i)
es.get_system<FEMSystem>(i).time = current_time;
}
void set_system_parameters(FEMSystem &system, FEMParameters ¶m)
{
system.verify_analytic_jacobians = param.verify_analytic_jacobians;
system.print_solution_norms = param.print_solution_norms;
system.print_solutions = param.print_solutions;
system.print_residual_norms = param.print_residual_norms;
system.print_residuals = param.print_residuals;
system.print_jacobian_norms = param.print_jacobian_norms;
system.print_jacobians = param.print_jacobians;
if (param.transient)
{
UnsteadySolver *innersolver;
if (param.timesolver_core == "euler2")
{
Euler2Solver *euler2solver =
new Euler2Solver(system);
euler2solver->theta = param.timesolver_theta;
innersolver = euler2solver;
}
else if (param.timesolver_core == "euler")
{
EulerSolver *eulersolver =
new EulerSolver(system);
eulersolver->theta = param.timesolver_theta;
innersolver = eulersolver;
}
else
{
std::cerr << "Don't recognize core TimeSolver type: "
<< param.timesolver_core << std::endl;
libmesh_error();
}
if (param.timesolver_tolerance)
{
TwostepTimeSolver *timesolver =
new TwostepTimeSolver(system);
timesolver->max_growth = param.timesolver_maxgrowth;
timesolver->target_tolerance = param.timesolver_tolerance;
timesolver->upper_tolerance = param.timesolver_upper_tolerance;
timesolver->component_norm = SystemNorm(param.timesolver_norm);
timesolver->core_time_solver =
AutoPtr<UnsteadySolver>(innersolver);
system.time_solver =
AutoPtr<UnsteadySolver>(timesolver);
}
else
system.time_solver =
AutoPtr<TimeSolver>(innersolver);
}
else
system.time_solver =
AutoPtr<TimeSolver>(new SteadySolver(system));
system.time_solver->reduce_deltat_on_diffsolver_failure =
param.deltat_reductions;
system.time_solver->quiet = param.time_solver_quiet;
system.deltat = param.deltat;
system.extra_quadrature_order = param.extra_quadrature_order;
if (param.use_petsc_snes)
{
#ifdef LIBMESH_HAVE_PETSC
PetscDiffSolver *solver = new PetscDiffSolver(system);
system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
#else
libmesh_error();
#endif
}
else
{
NewtonSolver *solver = new NewtonSolver(system);
system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
solver->quiet = param.solver_quiet;
solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
solver->minsteplength = param.min_step_length;
solver->relative_step_tolerance = param.relative_step_tolerance;
solver->relative_residual_tolerance = param.relative_residual_tolerance;
solver->require_residual_reduction = param.require_residual_reduction;
solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
if (system.time_solver->reduce_deltat_on_diffsolver_failure)
{
solver->continue_after_max_iterations = true;
solver->continue_after_backtrack_failure = true;
}
solver->max_linear_iterations = param.max_linear_iterations;
solver->initial_linear_tolerance = param.initial_linear_tolerance;
solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
}
}
#ifdef LIBMESH_ENABLE_AMR
AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
FEMParameters ¶m)
{
AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
mesh_refinement->coarsen_by_parents() = true;
mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
mesh_refinement->nelem_target() = param.nelem_target;
mesh_refinement->refine_fraction() = param.refine_fraction;
mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
return mesh_refinement;
}
#endif // LIBMESH_ENABLE_AMR
AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters& /* param */)
{
AutoPtr<ErrorEstimator> error_estimator;
error_estimator.reset(new KellyErrorEstimator);
return error_estimator;
}
AutoPtr<ErrorEstimator>
build_error_estimator_component_wise
(FEMParameters ¶m,
std::vector<std::vector<Real> > &term_weights,
std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type)
{
AutoPtr<ErrorEstimator> error_estimator;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
PatchRecoveryErrorEstimator *p1 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->primal_error_estimator().reset(p1);
PatchRecoveryErrorEstimator *p2 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->dual_error_estimator().reset(p2);
p1->set_patch_reuse(param.patch_reuse);
p2->set_patch_reuse(param.patch_reuse);
unsigned int size = primal_error_norm_type.size();
libmesh_assert_equal_to (size, dual_error_norm_type.size());
for (unsigned int i = 0; i != size; ++i)
{
adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
}
libmesh_assert_equal_to (size, term_weights.size());
for (unsigned int i = 0; i != term_weights.size(); ++i)
{
libmesh_assert_equal_to (size, term_weights[i].size());
adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
for (unsigned int j = 0; j != size; ++j)
if (i != j)
adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
}
return error_estimator;
}
AutoPtr<ErrorEstimator>
build_weighted_error_estimator_component_wise
(FEMParameters ¶m,
std::vector<std::vector<Real> > &term_weights,
std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type,
std::vector<FEMFunctionBase<Number>*> coupled_system_weight_functions)
{
AutoPtr<ErrorEstimator> error_estimator;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
WeightedPatchRecoveryErrorEstimator *p1 =
new WeightedPatchRecoveryErrorEstimator;
adjoint_residual_estimator->primal_error_estimator().reset(p1);
PatchRecoveryErrorEstimator *p2 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->dual_error_estimator().reset(p2);
p1->set_patch_reuse(param.patch_reuse);
p2->set_patch_reuse(param.patch_reuse);
p1->weight_functions.clear();
unsigned int size = primal_error_norm_type.size();
libmesh_assert(coupled_system_weight_functions.size() == size);
libmesh_assert(dual_error_norm_type.size() == size);
for (unsigned int i = 0; i != size; ++i)
{
p1->weight_functions.push_back(coupled_system_weight_functions[i]);
adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
}
libmesh_assert_equal_to (size, term_weights.size());
for (unsigned int i = 0; i != term_weights.size(); ++i)
{
libmesh_assert_equal_to (size, term_weights[i].size());
adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
for (unsigned int j = 0; j != size; ++j)
if (i != j)
adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
}
return error_estimator;
}
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
std::cout << "Started " << argv[0] << std::endl;
{
std::ifstream i("general.in");
if (!i)
{
std::cerr << '[' << libMesh::processor_id()
<< "] Can't find general.in; exiting early."
<< std::endl;
libmesh_error();
}
}
GetPot infile("general.in");
FEMParameters param;
param.read(infile);
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Mesh mesh (param.dimension);
AutoPtr<MeshRefinement> mesh_refinement =
build_mesh_refinement(mesh, param);
EquationSystems equation_systems (mesh);
std::cout << "Building mesh" << std::endl;
if (!param.initial_timestep && param.run_simulation)
build_domain(mesh, param);
std::cout << "Building system" << std::endl;
FEMSystem &system = build_system(equation_systems, infile, param);
set_system_parameters(system, param);
std::cout << "Initializing systems" << std::endl;
{
CoupledSystemQoI qoi;
qoi.assemble_qoi_sides = true;
system.attach_qoi(&qoi);
}
equation_systems.init ();
mesh.print_info();
equation_systems.print_info();
std::cout<<"Starting adaptive loop"<<std::endl<<std::endl;
unsigned int a_step = 0;
LinearSolver<Number> *linear_solver = system.get_linear_solver();
for (; a_step <= param.max_adaptivesteps; ++a_step)
{
if(param.global_tolerance > 0 && param.nelem_target > 0)
{
std::cout<<"We cant refine to both a non-zero tolerance and a target number of elements, EXITING adaptive loop. "<<std::endl<<std::endl;
break;
}
linear_solver->reuse_preconditioner(false);
std::cout<< "Adaptive step " << a_step << std::endl;
std::cout<< "Solving the forward problem" <<std::endl;
std::cout << "We have " << mesh.n_active_elem()
<< " active elements and " << equation_systems.n_active_dofs()
<< " active dofs." << std::endl << std::endl;
system.solve();
write_output(equation_systems, 0, a_step, "primal", param);
system.assemble_qoi();
system.postprocess();
Number QoI_0_computed = (dynamic_cast<CoupledSystem&>(system)).get_QoI_value();
std::cout<< "The boundary QoI is " << std::setprecision(17) << QoI_0_computed << std::endl << std::endl;
if(a_step!=param.max_adaptivesteps)
{
NumericVector<Number> &primal_solution = (*system.solution);
system.assemble_qoi_sides = true;
linear_solver->reuse_preconditioner(param.reuse_preconditioner);
std::cout<< "Solving the adjoint problem" <<std::endl;
system.adjoint_solve();
system.set_adjoint_already_solved(true);
NumericVector<Number> &dual_solution = system.get_adjoint_solution();
primal_solution.swap(dual_solution);
write_output(equation_systems, 0, a_step, "adjoint", param);
primal_solution.swap(dual_solution);
Real Pe = (dynamic_cast<CoupledSystem&>(system)).get_Pe();
ErrorVector error;
ErrorVector error_non_pressure;
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_non_pressure;
primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
primal_norm_type_vector_non_pressure.push_back(L2);
primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_non_pressure;
dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
dual_norm_type_vector_non_pressure.push_back(L2);
dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
std::vector<std::vector<Real> >
weights_matrix_non_pressure(system.n_vars(),
std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_non_pressure[0][0] = 1.;
weights_matrix_non_pressure[1][1] = 1.;
weights_matrix_non_pressure[3][3] = 1./Pe;
AutoPtr<ErrorEstimator> error_estimator_non_pressure =
build_error_estimator_component_wise
(param, weights_matrix_non_pressure,
primal_norm_type_vector_non_pressure,
dual_norm_type_vector_non_pressure);
error_estimator_non_pressure->estimate_error(system, error_non_pressure);
write_error(equation_systems, error_non_pressure, 0, a_step, param, "_non_pressure");
ErrorVector error_with_pressure;
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_with_pressure;
primal_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
primal_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
primal_norm_type_vector_with_pressure.push_back(L2);
primal_norm_type_vector_with_pressure.push_back(L2);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_with_pressure;
dual_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
dual_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
dual_norm_type_vector_with_pressure.push_back(L2);
dual_norm_type_vector_with_pressure.push_back(L2);
std::vector<std::vector<Real> >
weights_matrix_with_pressure
(system.n_vars(),
std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_with_pressure[0][2] = 1.;
weights_matrix_with_pressure[1][2] = 1.;
weights_matrix_with_pressure[2][0] = 1.;
weights_matrix_with_pressure[2][1] = 1.;
AutoPtr<ErrorEstimator> error_estimator_with_pressure =
build_error_estimator_component_wise
(param, weights_matrix_with_pressure,
primal_norm_type_vector_with_pressure,
dual_norm_type_vector_with_pressure);
error_estimator_with_pressure->estimate_error(system, error_with_pressure);
write_error(equation_systems, error_with_pressure, 0, a_step, param, "_with_pressure");
ErrorVector error_convection_diffusion_x;
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_convection_diffusion_x;
primal_norm_type_vector_convection_diffusion_x.push_back(L2);
primal_norm_type_vector_convection_diffusion_x.push_back(L2);
primal_norm_type_vector_convection_diffusion_x.push_back(L2);
primal_norm_type_vector_convection_diffusion_x.push_back(H1_X_SEMINORM);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_convection_diffusion_x;
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
dual_norm_type_vector_convection_diffusion_x.push_back(L2);
std::vector<std::vector<Real> >
weights_matrix_convection_diffusion_x
(system.n_vars(),
std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_convection_diffusion_x[0][3] = 1.;
weights_matrix_convection_diffusion_x[3][3] = 1.;
ConstFEMFunction<Number> identity(1);
CoupledFEMFunctionsx convdiffx0(system, 0);
CoupledFEMFunctionsx convdiffx3(system, 3);
std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_x;
coupled_system_weight_functions_x.push_back(&convdiffx0);
coupled_system_weight_functions_x.push_back(&identity);
coupled_system_weight_functions_x.push_back(&identity);
coupled_system_weight_functions_x.push_back(&convdiffx3);
AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_x =
build_weighted_error_estimator_component_wise
(param, weights_matrix_convection_diffusion_x,
primal_norm_type_vector_convection_diffusion_x,
dual_norm_type_vector_convection_diffusion_x,
coupled_system_weight_functions_x);
error_estimator_convection_diffusion_x->estimate_error
(system, error_convection_diffusion_x);
write_error(equation_systems, error_convection_diffusion_x,
0, a_step, param, "_convection_diffusion_x");
ErrorVector error_convection_diffusion_y;
std::vector<libMeshEnums::FEMNormType>
primal_norm_type_vector_convection_diffusion_y;
primal_norm_type_vector_convection_diffusion_y.push_back(L2);
primal_norm_type_vector_convection_diffusion_y.push_back(L2);
primal_norm_type_vector_convection_diffusion_y.push_back(L2);
primal_norm_type_vector_convection_diffusion_y.push_back(H1_Y_SEMINORM);
std::vector<libMeshEnums::FEMNormType>
dual_norm_type_vector_convection_diffusion_y;
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
dual_norm_type_vector_convection_diffusion_y.push_back(L2);
std::vector<std::vector<Real> >
weights_matrix_convection_diffusion_y
(system.n_vars(), std::vector<Real>(system.n_vars(), 0.0));
weights_matrix_convection_diffusion_y[1][3] = 1.;
weights_matrix_convection_diffusion_y[3][3] = 1.;
CoupledFEMFunctionsy convdiffy1(system, 1);
CoupledFEMFunctionsy convdiffy3(system, 3);
std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_y;
coupled_system_weight_functions_y.push_back(&identity);
coupled_system_weight_functions_y.push_back(&convdiffy1);
coupled_system_weight_functions_y.push_back(&identity);
coupled_system_weight_functions_y.push_back(&convdiffy3);
AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_y =
build_weighted_error_estimator_component_wise
(param, weights_matrix_convection_diffusion_y,
primal_norm_type_vector_convection_diffusion_y,
dual_norm_type_vector_convection_diffusion_y,
coupled_system_weight_functions_y);
error_estimator_convection_diffusion_y->estimate_error(system, error_convection_diffusion_y);
write_error(equation_systems, error_convection_diffusion_y, 0, a_step, param, "_convection_diffusion_y");
if(param.indicator_type == "adjoint_residual")
{
error.resize(error_non_pressure.size());
for(unsigned int i = 0; i < error.size(); i++)
{
error[i] = error_non_pressure[i] + error_with_pressure[i] + error_convection_diffusion_x[i] + error_convection_diffusion_y[i];
}
}
else
{
std::cout<<"Using Kelly Estimator"<<std::endl;
AutoPtr<ErrorEstimator> error_estimator = build_error_estimator(param);
error_estimator->estimate_error(system, error);
}
write_error(equation_systems, error, 0, a_step, param, "_total");
if(param.refine_uniformly)
{
mesh_refinement->uniformly_refine(1);
}
else if(param.global_tolerance >= 0. && param.nelem_target == 0.) // Refine based on reaching an error tolerance
{
mesh_refinement->flag_elements_by_error_tolerance (error);
mesh_refinement->refine_and_coarsen_elements();
}
else // Refine based on reaching a target number of elements
{
if (mesh.n_active_elem() >= param.nelem_target)
{
std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
break;
}
mesh_refinement->flag_elements_by_nelem_target (error);
mesh_refinement->refine_and_coarsen_elements();
}
equation_systems.reinit();
}
}
write_output_footers(param);
#endif // #ifndef LIBMESH_ENABLE_AMR
return 0;
}
The source file coupled_system.C without comments:
#include "libmesh/boundary_info.h"
#include "libmesh/dirichlet_boundaries.h"
#include "libmesh/dof_map.h"
#include "libmesh/fe_base.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fem_context.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh.h"
#include "libmesh/parallel.h"
#include "libmesh/quadrature.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/zero_function.h"
#include "coupled_system.h"
using namespace libMesh;
class BdyFunction : public FunctionBase<Number>
{
public:
BdyFunction (unsigned int u_var, unsigned int v_var, int sign)
: _u_var(u_var), _v_var(v_var), _sign(sign)
{ this->_initialized = true; }
virtual Number operator() (const Point&, const Real = 0)
{ libmesh_not_implemented(); }
virtual void operator() (const Point& p,
const Real,
DenseVector<Number>& output)
{
output.resize(2);
output.zero();
const Real y=p(1);
output(_u_var) = (_sign)*((y-2) * (y-3));
output(_v_var) = 0;
}
virtual AutoPtr<FunctionBase<Number> > clone() const
{ return AutoPtr<FunctionBase<Number> > (new BdyFunction(_u_var, _v_var, _sign)); }
private:
const unsigned int _u_var, _v_var;
const Real _sign;
};
void CoupledSystem::init_data ()
{
GetPot infile("coupled_system.in");
Peclet = infile("Peclet", 1.);
unsigned int pressure_p = infile("pressure_p", 1);
std::string fe_family = infile("fe_family", std::string("LAGRANGE"));
libmesh_assert((pressure_p == 1) || (fe_family != "LAGRANGE"));
FEFamily fefamily = Utility::string_to_enum<FEFamily>(fe_family);
u_var = this->add_variable ("u", static_cast<Order>(pressure_p+1),
fefamily);
v_var = this->add_variable ("v", static_cast<Order>(pressure_p+1),
fefamily);
p_var = this->add_variable ("p", static_cast<Order>(pressure_p),
fefamily);
C_var = this->add_variable ("C", static_cast<Order>(pressure_p+1),
fefamily);
this->time_evolving(u_var);
this->time_evolving(v_var);
this->time_evolving(C_var);
this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
this->print_jacobians = infile("print_jacobians", false);
this->print_element_jacobians = infile("print_element_jacobians", false);
const boundary_id_type left_inlet_id = 0;
std::set<boundary_id_type> left_inlet_bdy;
left_inlet_bdy.insert(left_inlet_id);
const boundary_id_type right_inlet_id = 1;
std::set<boundary_id_type> right_inlet_bdy;
right_inlet_bdy.insert(right_inlet_id);
const boundary_id_type outlets_id = 2;
std::set<boundary_id_type> outlets_bdy;
outlets_bdy.insert(outlets_id);
const boundary_id_type wall_id = 3;
std::set<boundary_id_type> wall_bdy;
wall_bdy.insert(wall_id);
std::vector<unsigned int> uv(1, u_var);
uv.push_back(v_var);
std::vector<unsigned int> C_only(1, C_var);
ZeroFunction<Number> zero;
ConstFunction<Number> one(1);
int velocity_sign = 1;
BdyFunction inflow_left(u_var, v_var, -velocity_sign);
BdyFunction inflow_right(u_var, v_var, velocity_sign);
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (wall_bdy, uv, &zero));
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (left_inlet_bdy, uv, &inflow_left));
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (left_inlet_bdy, C_only, &one));
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (right_inlet_bdy, uv, &inflow_right));
this->get_dof_map().add_dirichlet_boundary
(DirichletBoundary (right_inlet_bdy, C_only, &zero));
FEMSystem::init_data();
}
void CoupledSystem::init_context(DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
c.element_fe_var[u_var]->get_JxW();
c.element_fe_var[u_var]->get_phi();
c.element_fe_var[u_var]->get_dphi();
c.element_fe_var[u_var]->get_xyz();
c.element_fe_var[p_var]->get_phi();
c.side_fe_var[u_var]->get_JxW();
c.side_fe_var[u_var]->get_phi();
c.side_fe_var[u_var]->get_xyz();
}
bool CoupledSystem::element_time_derivative (bool request_jacobian,
DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW =
c.element_fe_var[u_var]->get_JxW();
const std::vector<std::vector<Real> >& phi =
c.element_fe_var[u_var]->get_phi();
const std::vector<std::vector<RealGradient> >& dphi =
c.element_fe_var[u_var]->get_dphi();
const std::vector<std::vector<Real> >& psi =
c.element_fe_var[p_var]->get_phi();
const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
libmesh_assert_equal_to (n_u_dofs, c.dof_indices_var[v_var].size());
DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
DenseSubMatrix<Number> &Kup = *c.elem_subjacobians[u_var][p_var];
DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
DenseSubMatrix<Number> &Kvv = *c.elem_subjacobians[v_var][v_var];
DenseSubMatrix<Number> &Kvp = *c.elem_subjacobians[v_var][p_var];
DenseSubVector<Number> &Fv = *c.elem_subresiduals[v_var];
DenseSubMatrix<Number> &KCu = *c.elem_subjacobians[C_var][u_var];
DenseSubMatrix<Number> &KCv = *c.elem_subjacobians[C_var][v_var];
DenseSubMatrix<Number> &KCC = *c.elem_subjacobians[C_var][C_var];
DenseSubVector<Number> &FC = *c.elem_subresiduals[C_var];
unsigned int n_qpoints = c.element_qrule->n_points();
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Number p = c.interior_value(p_var, qp),
u = c.interior_value(u_var, qp),
v = c.interior_value(v_var, qp);
Gradient grad_u = c.interior_gradient(u_var, qp),
grad_v = c.interior_gradient(v_var, qp),
grad_C = c.interior_gradient(C_var, qp);
NumberVectorValue U (u, v);
const Number C_x = grad_C(0);
const Number C_y = grad_C(1);
for (unsigned int i=0; i != n_u_dofs; i++)
{
Fu(i) += JxW[qp] *
(p*dphi[i][qp](0) - // pressure term
(grad_u*dphi[i][qp])); // diffusion term
Fv(i) += JxW[qp] *
(p*dphi[i][qp](1) - // pressure term
(grad_v*dphi[i][qp])); // diffusion term
FC(i) += JxW[qp] *
( (U*grad_C)*phi[i][qp] + // convection term
(1./Peclet)*(grad_C*dphi[i][qp]) ); // diffusion term
if (request_jacobian && c.elem_solution_derivative)
{
libmesh_assert (c.elem_solution_derivative == 1.0);
for (unsigned int j=0; j != n_u_dofs; j++)
{
Kuu(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term */
Kvv(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term */
KCu(i,j) += JxW[qp]* ( (phi[j][qp]*C_x)*phi[i][qp] ); /* convection term */
KCv(i,j) += JxW[qp]*( (phi[j][qp]*C_y)*phi[i][qp] ); /* convection term */
KCC(i,j) += JxW[qp]*
( (U*dphi[j][qp])*phi[i][qp] + /* nonlinear term (convection) */
(1./Peclet)*(dphi[j][qp]*dphi[i][qp]) ); /* diffusion term */
}
for (unsigned int j=0; j != n_p_dofs; j++)
{
Kup(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](0);
Kvp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](1);
}
}
}
} // end of the quadrature point qp-loop
return request_jacobian;
}
bool CoupledSystem::element_constraint (bool request_jacobian,
DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.element_fe_var[u_var]->get_JxW();
const std::vector<std::vector<RealGradient> >& dphi =
c.element_fe_var[u_var]->get_dphi();
const std::vector<std::vector<Real> >& psi =
c.element_fe_var[p_var]->get_phi();
const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
DenseSubMatrix<Number> &Kpu = *c.elem_subjacobians[p_var][u_var];
DenseSubMatrix<Number> &Kpv = *c.elem_subjacobians[p_var][v_var];
DenseSubVector<Number> &Fp = *c.elem_subresiduals[p_var];
unsigned int n_qpoints = c.element_qrule->n_points();
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Gradient grad_u = c.interior_gradient(u_var, qp),
grad_v = c.interior_gradient(v_var, qp);
for (unsigned int i=0; i != n_p_dofs; i++)
{
Fp(i) += JxW[qp] * psi[i][qp] *
(grad_u(0) + grad_v(1));
if (request_jacobian && c.elem_solution_derivative)
{
libmesh_assert (c.elem_solution_derivative == 1.0);
for (unsigned int j=0; j != n_u_dofs; j++)
{
Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
}
}
}
} // end of the quadrature point qp-loop
return request_jacobian;
}
void CoupledSystem::postprocess()
{
computed_QoI = 0.0;
computed_QoI = System::qoi[0];
}
Number CoupledFEMFunctionsx::operator()(const FEMContext& c, const Point& p,
const Real /* time */)
{
Number weight = 0.0;
switch(var)
{
case 0:
{
Gradient grad_C = c.point_gradient(3, p);
weight = grad_C(0);
}
break;
case 3:
{
Number u = c.point_value(0, p);
weight = u;
}
break;
default:
{
std::cout<<"Wrong variable number"<<var<<" passed to CoupledFEMFunctionsx object ! Quitting !"<<std::endl;
libmesh_error();
}
}
return weight;
}
Number CoupledFEMFunctionsy::operator()(const FEMContext& c, const Point& p,
const Real /* time */)
{
Number weight = 0.0;
switch(var)
{
case 1:
{
Gradient grad_C = c.point_gradient(3, p);
weight = grad_C(1);
}
break;
case 3:
{
Number v = c.point_value(1, p);
weight = v;
}
break;
default:
{
std::cout<<"Wrong variable number "<<var<<" passed to CoupledFEMFunctionsy object ! Quitting !"<<std::endl;
libmesh_error();
}
}
return weight;
}
The source file domain.C without comments:
#include "libmesh/mesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_refinement.h"
#include "domain.h"
using namespace libMesh;
void build_domain (Mesh &mesh, FEMParameters ¶m)
{
mesh.read(param.domainfile);
std::cout<<"Making elements 2nd order"<<std::endl;
mesh.all_second_order();
}
The source file femparameters.C without comments:
#include "femparameters.h"
using namespace libMesh;
#define GETPOT_INPUT(A) { A = input(#A, A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_INT_INPUT(A) { A = input(#A, (int)A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_REGISTER(A) { \
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
void FEMParameters::read(GetPot &input)
{
std::vector<std::string> variable_names;
GETPOT_INT_INPUT(initial_timestep);
GETPOT_INT_INPUT(n_timesteps);
GETPOT_INPUT(transient);
GETPOT_INT_INPUT(deltat_reductions);
GETPOT_INPUT(timesolver_core);
GETPOT_INPUT(end_time);
GETPOT_INPUT(deltat);
GETPOT_INPUT(timesolver_theta);
GETPOT_INPUT(timesolver_maxgrowth);
GETPOT_INPUT(timesolver_upper_tolerance);
GETPOT_INPUT(timesolver_tolerance);
GETPOT_INPUT(steadystate_tolerance);
GETPOT_REGISTER(timesolver_norm);
const unsigned int n_timesolver_norm = input.vector_variable_size("timesolver_norm");
timesolver_norm.resize(n_timesolver_norm, L2);
for (unsigned int i=0; i != n_timesolver_norm; ++i)
{
int current_norm = 0; // L2
if (timesolver_norm[i] == H1)
current_norm = 1;
if (timesolver_norm[i] == H2)
current_norm = 2;
current_norm = input("timesolver_norm", current_norm, i);
if (current_norm == 0)
timesolver_norm[i] = L2;
else if (current_norm == 1)
timesolver_norm[i] = H1;
else if (current_norm == 2)
timesolver_norm[i] = H2;
else
timesolver_norm[i] = DISCRETE_L2;
}
GETPOT_INT_INPUT(dimension);
GETPOT_INPUT(domaintype);
GETPOT_INPUT(domainfile);
GETPOT_INPUT(elementtype);
GETPOT_INPUT(fine_mesh_file_primal);
GETPOT_INPUT(fine_mesh_soln_primal);
GETPOT_INPUT(fine_mesh_file_adjoint);
GETPOT_INPUT(fine_mesh_soln_adjoint);
GETPOT_INPUT(elementorder);
GETPOT_INPUT(domain_xmin);
GETPOT_INPUT(domain_ymin);
GETPOT_INPUT(domain_zmin);
GETPOT_INPUT(domain_edge_width);
GETPOT_INPUT(domain_edge_length);
GETPOT_INPUT(domain_edge_height);
GETPOT_INT_INPUT(coarsegridx);
GETPOT_INT_INPUT(coarsegridy);
GETPOT_INT_INPUT(coarsegridz);
GETPOT_INT_INPUT(coarserefinements);
GETPOT_INT_INPUT(coarsecoarsenings);
GETPOT_INT_INPUT(extrarefinements);
GETPOT_INPUT(use_petsc_snes);
GETPOT_INPUT(time_solver_quiet);
GETPOT_INPUT(solver_quiet);
GETPOT_INPUT(reuse_preconditioner);
GETPOT_INPUT(require_residual_reduction);
GETPOT_INPUT(min_step_length);
GETPOT_INT_INPUT(max_linear_iterations);
GETPOT_INT_INPUT(max_nonlinear_iterations);
GETPOT_INPUT(relative_step_tolerance);
GETPOT_INPUT(relative_residual_tolerance);
GETPOT_INPUT(initial_linear_tolerance);
GETPOT_INPUT(minimum_linear_tolerance);
GETPOT_INPUT(linear_tolerance_multiplier);
GETPOT_INT_INPUT(nelem_target);
GETPOT_INPUT(global_tolerance);
GETPOT_INPUT(refine_fraction);
GETPOT_INPUT(coarsen_fraction);
GETPOT_INPUT(coarsen_threshold);
GETPOT_INT_INPUT(max_adaptivesteps);
GETPOT_INT_INPUT(initial_adaptivesteps);
GETPOT_INT_INPUT(initial_sobolev_order);
GETPOT_INT_INPUT(initial_extra_quadrature);
GETPOT_INPUT(refine_uniformly);
GETPOT_INPUT(indicator_type);
GETPOT_INPUT(adjoint_residual_type);
GETPOT_INPUT(patch_reuse);
GETPOT_INT_INPUT(sobolev_order);
GETPOT_INPUT(alternate_with_uniform_steps);
GETPOT_INPUT(alternate_step_number);
GETPOT_INPUT(component_wise_error);
GETPOT_INPUT(compute_sensitivities);
GETPOT_INPUT(compare_to_fine_solution_primal);
GETPOT_INPUT(compare_to_fine_solution_adjoint);
GETPOT_INPUT(do_forward_sensitivity);
GETPOT_INPUT(do_adjoint_sensitivity);
GETPOT_INPUT(postprocess_adjoint);
GETPOT_INT_INPUT(write_interval);
GETPOT_INPUT(output_xda);
GETPOT_INPUT(output_xdr);
GETPOT_INPUT(output_gz);
#ifndef LIBMESH_HAVE_GZSTREAM
output_gz = false;
#endif
GETPOT_INPUT(output_bz2);
#ifndef LIBMESH_HAVE_BZ2
output_bz2 = false;
#endif
GETPOT_INPUT(output_gmv);
GETPOT_INPUT(write_gmv_error);
#ifndef LIBMESH_HAVE_GMV
output_gmv = false;
write_gmv_error = false;
#endif
GETPOT_INPUT(output_tecplot);
GETPOT_INPUT(write_tecplot_error);
#ifndef LIBMESH_HAVE_TECPLOT_API
output_tecplot = false;
write_tecplot_error = false;
#endif
GETPOT_INPUT(run_simulation);
GETPOT_INPUT(run_postprocess);
run_simulation = input("run_simulation", run_simulation);
run_postprocess = input("run_postprocess", run_postprocess);
GETPOT_REGISTER(fe_family);
const unsigned int n_fe_family =
std::max(1u, input.vector_variable_size("fe_family"));
fe_family.resize(n_fe_family, "LAGRANGE");
for (unsigned int i=0; i != n_fe_family; ++i)
fe_family[i] = input("fe_family", fe_family[i].c_str(), i);
GETPOT_REGISTER(fe_order);
const unsigned int n_fe_order =
input.vector_variable_size("fe_order");
fe_order.resize(n_fe_order, 1);
for (unsigned int i=0; i != n_fe_order; ++i)
fe_order[i] = input("fe_order", (int)fe_order[i], i);
GETPOT_INPUT(extra_quadrature_order);
GETPOT_INPUT(analytic_jacobians);
GETPOT_INPUT(verify_analytic_jacobians);
GETPOT_INPUT(print_solution_norms);
GETPOT_INPUT(print_solutions);
GETPOT_INPUT(print_residual_norms);
GETPOT_INPUT(print_residuals);
GETPOT_INPUT(print_jacobian_norms);
GETPOT_INPUT(print_jacobians);
}
The source file H-qoi.C without comments:
#include "H-qoi.h"
using namespace libMesh;
void CoupledSystemQoI::init_qoi( std::vector<Number>& sys_qoi)
{
sys_qoi.resize(1);
return;
}
void CoupledSystemQoI::side_qoi_derivative (DiffContext &context,
const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
const std::vector<std::vector<Real> > &phi = c.side_fe_var[0]->get_phi();
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
const unsigned int n_u_dofs = c.dof_indices_var[1].size();
DenseSubVector<Number> &Qu = *c.elem_qoi_subderivatives[0][0];
DenseSubVector<Number> &QC = *c.elem_qoi_subderivatives[0][3];
unsigned int n_qpoints = c.side_qrule->n_points();
Number u = 0. ;
Number C = 0. ;
if(c.has_side_boundary_id(2))
{
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Real x = q_point[qp](0);
if(x < 0.)
{
u = c.side_value(0,qp);
C = c.side_value(3,qp);
for (unsigned int i=0; i != n_u_dofs; i++)
{
Qu(i) += JxW[qp] * -phi[i][qp] * C;
QC(i) += JxW[qp] * phi[i][qp] * -u;
}
} // end if
} // end quadrature loop
} // end if on outlet
}
void CoupledSystemQoI::side_qoi(DiffContext &context, const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
unsigned int n_qpoints = c.side_qrule->n_points();
Number dQoI_0 = 0. ;
Number u = 0. ;
Number C = 0. ;
if(c.has_side_boundary_id(2))
{
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Real x = q_point[qp](0);
if(x < 0.)
{
u = c.side_value(0,qp);
C = c.side_value(3,qp);
dQoI_0 += JxW[qp] * -u * C;
} // end if
} // end quadrature loop
} // end if on bdry
c.elem_qoi[0] += dQoI_0;
}
The source file initial.C without comments:
#include "initial.h"
#include "femparameters.h"
using namespace libMesh;
void read_initial_parameters()
{
}
void finish_initialization()
{
}
Number initial_value(const Point& /* p */,
const Parameters& /* param */,
const std::string&,
const std::string&)
{
return 1.;
}
Gradient initial_grad(const Point& /* p */,
const Parameters& /* param */,
const std::string&,
const std::string&)
{
return 0.;
}
The source file output.C without comments:
#include <fstream>
#include <unistd.h>
#include "libmesh/libmesh_common.h"
#include "output.h"
using namespace libMesh;
void start_output(unsigned int timesteps,
std::string filename,
std::string varname)
{
if (timesteps)
{
std::ifstream infile (filename.c_str());
char buf[1025];
if (!infile.is_open())
libmesh_error();
for (unsigned int i=0; i != timesteps+1; ++i)
infile.getline(buf, 1024);
if (!infile.good())
libmesh_error();
unsigned int length = infile.tellg();
infile.close();
int err = truncate(filename.c_str(), length);
if (err != 0)
libmesh_error();
}
else
{
std::ofstream outfile (filename.c_str(), std::ios_base::out);
outfile << varname << " = [" << std::endl;
}
}
The console output of the program:
***************************************************************
* Running Example adjoints_ex3:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Started /workspace/libmesh/examples/adjoints/adjoints_ex3/.libs/lt-example-devel
Building mesh
Making elements 2nd order
Building system
Initializing systems
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=99
n_local_nodes()=51
n_elem()=16
n_local_elem()=8
n_active_elem()=16
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "CoupledSystem"
Type "Implicit"
Variables={ "u" "v" } "p" "C"
Finite Element Types="LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN" "CARTESIAN" "CARTESIAN"
Approximation Orders="SECOND", "THIRD" "FIRST", "THIRD" "SECOND", "THIRD"
n_dofs()=331
n_local_dofs()=171
n_constrained_dofs()=138
n_local_constrained_dofs()=66
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 40.006
Average Off-Processor Bandwidth <= 1.32931
Maximum On-Processor Bandwidth <= 71
Maximum Off-Processor Bandwidth <= 20
DofMap Constraints
Number of DoF Constraints = 138
Number of Heterogenous Constraints= 5
Average DoF Constraint Length= 0
Number of Node Constraints = 0
Starting adaptive loop
Adaptive step 0
Solving the forward problem
We have 16 active elements and 193 active dofs.
Nonlinear solver converged, step 1, residual reduction 7.65258e-14 < 1e-05
The boundary QoI is 0.083342124064179859
Solving the adjoint problem
Adaptive step 1
Solving the forward problem
We have 22 active elements and 259 active dofs.
Nonlinear solver converged, step 0, residual reduction 5.6868336442060613e-07 < 1.0000000000000001e-05
The boundary QoI is 0.08228698767295109
Solving the adjoint problem
Adaptive step 2
Solving the forward problem
We have 31 active elements and 363 active dofs.
Nonlinear solver converged, step 1, residual reduction 4.0807919596899815e-14 < 1.0000000000000001e-05
Nonlinear solver converged, step 1, relative step size 7.5013669085787913e-07 < 1.0000000000000001e-05
The boundary QoI is 0.083169669189009976
Solving the adjoint problem
Adaptive step 3
Solving the forward problem
We have 43 active elements and 509 active dofs.
Nonlinear solver converged, step 1, residual reduction 3.4075021297015425e-14 < 1.0000000000000001e-05
Nonlinear solver converged, step 1, relative step size 6.4968541447564801e-06 < 1.0000000000000001e-05
The boundary QoI is 0.083351743975040624
Solving the adjoint problem
Adaptive step 4
Solving the forward problem
We have 58 active elements and 681 active dofs.
Nonlinear solver converged, step 0, residual reduction 5.1805208687554768e-06 < 1.0000000000000001e-05
The boundary QoI is 0.08333729860703705
Solving the adjoint problem
Adaptive step 5
Solving the forward problem
We have 76 active elements and 908 active dofs.
Nonlinear solver converged, step 0, residual reduction 6.4317035422377305e-07 < 1.0000000000000001e-05
The boundary QoI is 0.083251251683180455
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adjoints/adjoints_ex3/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:39:51 2013
Using Petsc Release Version 3.3.0, Patch 2, Fri Jul 13 15:42:00 CDT 2012
Max Max/Min Avg Total
Time (sec): 1.234e+01 1.00000 1.234e+01
Objects: 8.980e+02 1.00000 8.980e+02
Flops: 1.355e+08 1.25950 1.215e+08 2.430e+08
Flops/sec: 1.098e+07 1.25950 9.852e+06 1.970e+07
MPI Messages: 7.320e+02 1.00000 7.320e+02 1.464e+03
MPI Message Lengths: 3.851e+05 1.00000 5.260e+02 7.701e+05
MPI Reductions: 2.468e+03 1.00000
Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N --> 2N flops
and VecAXPY() for complex vectors of length N --> 8N flops
Summary of Stages: ----- Time ------ ----- Flops ----- --- Messages --- -- Message Lengths -- -- Reductions --
Avg %Total Avg %Total counts %Total Avg %Total counts %Total
0: Main Stage: 1.2335e+01 100.0% 2.4305e+08 100.0% 1.464e+03 100.0% 5.260e+02 100.0% 2.467e+03 100.0%
------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
Count: number of times phase was executed
Time and Flops: Max - maximum over all processors
Ratio - ratio of maximum to minimum over all processors
Mess: number of messages sent
Avg. len: average message length
Reduct: number of global reductions
Global: entire computation
Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
%T - percent time in this phase %f - percent flops in this phase
%M - percent messages in this phase %L - percent message lengths in this phase
%R - percent reductions in this phase
Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event Count Time (sec) Flops --- Global --- --- Stage --- Total
Max Ratio Max Ratio Max Ratio Mess Avg len Reduct %T %f %M %L %R %T %f %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------
--- Event Stage 0: Main Stage
KSPGMRESOrthog 235 1.0 4.5779e-03 2.0 3.50e+06 1.1 0.0e+00 0.0e+00 2.4e+02 0 3 0 0 10 0 3 0 0 10 1490
KSPSetUp 28 1.0 2.1505e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
KSPSolve 9 1.0 2.2387e-01 1.0 1.04e+08 1.3 3.8e+02 2.4e+02 4.5e+02 2 75 26 12 18 2 75 26 12 18 817
PCSetUp 28 1.0 2.6817e-01 1.2 8.06e+07 1.4 0.0e+00 0.0e+00 1.2e+02 2 57 0 0 5 2 57 0 0 5 520
PCSetUpOnBlocks 9 1.0 1.9072e-01 1.3 5.87e+07 1.5 0.0e+00 0.0e+00 6.3e+01 1 41 0 0 3 1 41 0 0 3 517
PCApply 264 1.0 9.3940e-02 1.0 6.41e+07 1.1 0.0e+00 0.0e+00 3.5e+01 1 49 0 0 1 1 49 0 0 1 1279
VecMDot 235 1.0 3.8898e-03 2.5 1.75e+06 1.1 0.0e+00 0.0e+00 2.4e+02 0 1 0 0 10 0 1 0 0 10 876
VecNorm 300 1.0 7.4594e-03 1.3 2.14e+05 1.1 0.0e+00 0.0e+00 3.0e+02 0 0 0 0 12 0 0 0 0 12 56
VecScale 250 1.0 1.2326e-04 1.0 9.01e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1424
VecCopy 244 1.0 1.4734e-04 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecSet 628 1.0 2.1529e-04 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecAXPY 39 1.0 5.3406e-05 1.1 2.75e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1001
VecMAXPY 250 1.0 5.7125e-04 1.0 1.92e+06 1.1 0.0e+00 0.0e+00 0.0e+00 0 2 0 0 0 0 2 0 0 0 6551
VecAssemblyBegin 409 1.0 6.1320e-02 1.8 0.00e+00 0.0 6.6e+01 1.2e+02 1.1e+03 0 0 5 1 43 0 0 5 1 43 0
VecAssemblyEnd 409 1.0 2.0146e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecScatterBegin 515 1.0 8.3518e-04 1.0 0.00e+00 0.0 1.0e+03 4.8e+02 0.0e+00 0 0 70 64 0 0 0 70 64 0 0
VecScatterEnd 515 1.0 5.3260e-02 4.4 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecNormalize 250 1.0 7.4885e-03 1.2 2.70e+05 1.1 0.0e+00 0.0e+00 2.5e+02 0 0 0 0 10 0 0 0 0 10 70
MatMult 188 1.0 5.6099e-02 3.7 6.85e+06 1.1 3.8e+02 2.4e+02 0.0e+00 0 5 26 12 0 0 5 26 12 0 235
MatMultTranspose 62 1.0 1.3418e-03 1.1 1.78e+06 1.1 1.2e+02 2.0e+02 0.0e+00 0 1 8 3 0 0 1 8 3 0 2531
MatSolve 197 1.0 1.2407e-02 1.2 3.47e+07 1.2 0.0e+00 0.0e+00 0.0e+00 0 27 0 0 0 0 27 0 0 0 5193
MatSolveTranspos 67 1.0 3.7863e-03 1.1 7.59e+06 1.1 0.0e+00 0.0e+00 0.0e+00 0 6 0 0 0 0 6 0 0 0 3884
MatLUFactorNum 14 1.0 4.1730e-02 1.3 8.06e+07 1.4 0.0e+00 0.0e+00 0.0e+00 0 57 0 0 0 0 57 0 0 0 3344
MatILUFactorSym 14 1.0 2.2373e-01 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 4.2e+01 2 0 0 0 2 2 0 0 0 2 0
MatAssemblyBegin 28 1.0 8.9536e-03 2.0 0.00e+00 0.0 4.2e+01 3.6e+03 5.6e+01 0 0 3 19 2 0 0 3 19 2 0
MatAssemblyEnd 28 1.0 1.5202e-03 1.2 0.00e+00 0.0 2.4e+01 6.0e+01 4.8e+01 0 0 2 0 2 0 0 2 0 2 0
MatGetRowIJ 14 1.0 4.7684e-06 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatGetOrdering 14 1.0 6.2251e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 5.6e+01 0 0 0 0 2 0 0 0 0 2 0
MatZeroEntries 26 1.0 2.6393e-04 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
------------------------------------------------------------------------------------------------------------------------
Memory usage is given in bytes:
Object Type Creations Destructions Memory Descendants' Mem.
Reports information only for process 0.
--- Event Stage 0: Main Stage
Krylov Solver 23 23 214192 0
Preconditioner 23 23 20448 0
Vector 515 515 2019784 0
Vector Scatter 83 83 85988 0
Index Set 190 190 190808 0
IS L to G Mapping 31 31 17484 0
Matrix 32 32 13130632 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 2.00272e-06
Average time for zero size MPI_Send(): 1.29938e-05
#PETSc Option Table entries:
-ksp_right_pc
-log_summary
-pc_type bjacobi
-sub_pc_factor_levels 4
-sub_pc_factor_zeropivot 0
-sub_pc_type ilu
#End of PETSc Option Table entries
Compiled without FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 sizeof(PetscInt) 4
Configure run at: Thu Nov 8 11:21:02 2012
Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared-libraries=1 --with-mpi-dir=/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1 --with-mumps=true --download-mumps=1 --with-metis=true --download-metis=1 --with-parmetis=true --download-parmetis=1 --with-superlu=true --download-superlu=1 --with-superludir=true --download-superlu_dist=1 --with-blacs=true --download-blacs=1 --with-scalapack=true --download-scalapack=1 --with-hypre=true --download-hypre=1 --with-blas-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_intel_lp64.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_sequential.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_core.so]" --with-lapack-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_lapack95_lp64.a]"
-----------------------------------------
Libraries compiled on Thu Nov 8 11:21:02 2012 on daedalus.ices.utexas.edu
Machine characteristics: Linux-2.6.32-279.1.1.el6.x86_64-x86_64-with-redhat-6.3-Carbon
Using PETSc directory: /opt/apps/ossw/libraries/petsc/petsc-3.3-p2
Using PETSc arch: intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt
-----------------------------------------
Using C compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx -wd1572 -O3 -fPIC ${COPTFLAGS} ${CFLAGS}
Using Fortran compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90 -fPIC -O3 ${FOPTFLAGS} ${FFLAGS}
-----------------------------------------
Using include paths: -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/include
-----------------------------------------
Using C linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx
Using Fortran linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90
Using libraries: -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lpetsc -lX11 -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lHYPRE -lpthread -lsuperlu_dist_3.0 -lparmetis -lmetis -lscalapack -lblacs -lsuperlu_4.3 -Wl,-rpath,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -L/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,-rpath,/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -L/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -Wl,-rpath,/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -L/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -L/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -lmpichf90 -lifport -lifcore -lm -lm -lmpichcxx -ldl -lmpich -lopa -lmpl -lrt -lpthread -limf -lsvml -lipgo -ldecimal -lcilkrts -lstdc++ -lgcc_s -lirc -lirc_s -ldl
-----------------------------------------
----------------------------------------------------------------------------------------------------------------------
| Processor id: 0 |
| Num Processors: 2 |
| Time: Tue Feb 5 13:39:51 2013 |
| OS: Linux |
| HostName: hbar.ices.utexas.edu |
| OS Release: 2.6.32-279.1.1.el6.x86_64 |
| OS Version: #1 SMP Tue Jul 10 11:24:23 CDT 2012 |
| Machine: x86_64 |
| Username: benkirk |
| Configuration: ./configure '--enable-everything' |
| '--prefix=/workspace/libmesh/install' |
| 'CXX=mpicxx' |
| 'CC=mpicc' |
| 'F77=mpif77' |
| 'FC=mpif90' |
| 'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2' |
| 'PETSC_ARCH=intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt' |
| 'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt' |
| 'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/intel-12.1/mpich2-1.4.1p1/mkl-intel-10.3.12.361'|
| 'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/intel-12.1' |
----------------------------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=12.3466, Active time=12.2328 |
---------------------------------------------------------------------------------------------------------------------
| Event nCalls Total Time Avg Time Total Time Avg Time % of Active Time |
| w/o Sub w/o Sub With Sub With Sub w/o S With S |
|---------------------------------------------------------------------------------------------------------------------|
| |
| |
| AdjointResidualErrorEstimator |
| estimate_error() 20 0.0059 0.000293 9.2526 0.462628 0.05 75.64 |
| |
| DofMap |
| add_neighbors_to_send_list() 31 0.0405 0.001306 0.0503 0.001624 0.33 0.41 |
| build_constraint_matrix() 516 0.0267 0.000052 0.0267 0.000052 0.22 0.22 |
| build_sparsity() 6 0.0522 0.008706 0.0978 0.016292 0.43 0.80 |
| cnstrn_elem_mat_vec() 171 0.0167 0.000098 0.0167 0.000098 0.14 0.14 |
| constrain_elem_matrix() 87 0.0035 0.000041 0.0035 0.000041 0.03 0.03 |
| constrain_elem_vector() 258 0.0038 0.000015 0.0038 0.000015 0.03 0.03 |
| create_dof_constraints() 31 0.0686 0.002214 0.3839 0.012383 0.56 3.14 |
| distribute_dofs() 31 0.0359 0.001157 0.0858 0.002767 0.29 0.70 |
| dof_indices() 64989 5.7115 0.000088 5.7115 0.000088 46.69 46.69 |
| enforce_constraints_exactly() 36 0.0062 0.000172 0.0062 0.000172 0.05 0.05 |
| old_dof_indices() 1180 0.1549 0.000131 0.1549 0.000131 1.27 1.27 |
| prepare_send_list() 31 0.0003 0.000009 0.0003 0.000009 0.00 0.00 |
| reinit() 31 0.0443 0.001429 0.0443 0.001429 0.36 0.36 |
| |
| EquationSystems |
| build_discontinuous_solution_vector() 25 0.0028 0.000111 0.0116 0.000462 0.02 0.09 |
| build_solution_vector() 11 0.0049 0.000442 0.0812 0.007381 0.04 0.66 |
| |
| FE |
| compute_shape_functions() 61100 1.1709 0.000019 1.1709 0.000019 9.57 9.57 |
| init_shape_functions() 9612 0.3873 0.000040 0.3873 0.000040 3.17 3.17 |
| inverse_map() 19356 0.1753 0.000009 0.1753 0.000009 1.43 1.43 |
| |
| FEMSystem |
| assemble_qoi() 6 0.0194 0.003236 0.2692 0.044874 0.16 2.20 |
| assemble_qoi_derivative() 5 0.0208 0.004169 0.2090 0.041798 0.17 1.71 |
| assembly() 9 0.1484 0.016492 0.4748 0.052751 1.21 3.88 |
| assembly(get_jacobian) 5 0.0744 0.014885 0.2390 0.047798 0.61 1.95 |
| assembly(get_residual) 9 0.0473 0.005259 0.3584 0.039817 0.39 2.93 |
| |
| FEMap |
| compute_affine_map() 61100 0.5871 0.000010 0.5871 0.000010 4.80 4.80 |
| compute_face_map() 3560 0.0822 0.000023 0.2140 0.000060 0.67 1.75 |
| init_face_shape_functions() 804 0.0040 0.000005 0.0040 0.000005 0.03 0.03 |
| init_reference_to_physical_map() 9612 0.3901 0.000041 0.3901 0.000041 3.19 3.19 |
| |
| GMVIO |
| write_nodal_data() 11 0.0110 0.001002 0.0116 0.001054 0.09 0.09 |
| |
| ImplicitSystem |
| adjoint_solve() 5 0.0004 0.000089 0.5426 0.108523 0.00 4.44 |
| |
| LocationMap |
| find() 400 0.0017 0.000004 0.0017 0.000004 0.01 0.01 |
| init() 10 0.0021 0.000211 0.0021 0.000211 0.02 0.02 |
| |
| Mesh |
| all_first_order() 25 0.0116 0.000464 0.0116 0.000464 0.09 0.09 |
| all_second_order() 1 0.0012 0.001171 0.0012 0.001171 0.01 0.01 |
| contract() 5 0.0001 0.000017 0.0002 0.000044 0.00 0.00 |
| find_neighbors() 57 0.0245 0.000431 0.0277 0.000485 0.20 0.23 |
| renumber_nodes_and_elem() 69 0.0014 0.000020 0.0014 0.000020 0.01 0.01 |
| |
| MeshCommunication |
| compute_hilbert_indices() 32 0.0054 0.000168 0.0054 0.000168 0.04 0.04 |
| find_global_indices() 32 0.0036 0.000112 0.0149 0.000466 0.03 0.12 |
| parallel_sort() 32 0.0032 0.000101 0.0041 0.000127 0.03 0.03 |
| |
| MeshOutput |
| write_equation_systems() 11 0.0004 0.000035 0.0937 0.008522 0.00 0.77 |
| |
| MeshRefinement |
| _coarsen_elements() 10 0.0001 0.000014 0.0002 0.000017 0.00 0.00 |
| _refine_elements() 10 0.0022 0.000220 0.0055 0.000555 0.02 0.05 |
| add_point() 400 0.0015 0.000004 0.0032 0.000008 0.01 0.03 |
| make_coarsening_compatible() 20 0.0015 0.000074 0.0015 0.000074 0.01 0.01 |
| make_refinement_compatible() 20 0.0002 0.000010 0.0003 0.000013 0.00 0.00 |
| |
| MetisPartitioner |
| partition() 32 0.0158 0.000495 0.0334 0.001043 0.13 0.27 |
| |
| NewtonSolver |
| solve() 6 0.0254 0.004233 1.0909 0.181811 0.21 8.92 |
| |
| Parallel |
| allgather() 169 0.0014 0.000008 0.0017 0.000010 0.01 0.01 |
| broadcast() 22 0.0004 0.000016 0.0003 0.000013 0.00 0.00 |
| max(bool) 52 0.0002 0.000003 0.0002 0.000003 0.00 0.00 |
| max(scalar) 3936 0.0105 0.000003 0.0105 0.000003 0.09 0.09 |
| max(vector) 900 0.0056 0.000006 0.0127 0.000014 0.05 0.10 |
| max(vector) 10 0.0001 0.000014 0.0002 0.000024 0.00 0.00 |
| min(bool) 4727 0.0117 0.000002 0.0117 0.000002 0.10 0.10 |
| min(scalar) 3819 0.2083 0.000055 0.2083 0.000055 1.70 1.70 |
| min(vector) 905 0.0063 0.000007 0.0138 0.000015 0.05 0.11 |
| probe() 316 0.0010 0.000003 0.0010 0.000003 0.01 0.01 |
| receive() 316 0.0013 0.000004 0.0024 0.000007 0.01 0.02 |
| send() 316 0.0007 0.000002 0.0007 0.000002 0.01 0.01 |
| send_receive() 380 0.0018 0.000005 0.0055 0.000014 0.01 0.04 |
| sum() 353 0.0025 0.000007 0.0138 0.000039 0.02 0.11 |
| |
| Parallel::Request |
| wait() 316 0.0005 0.000001 0.0005 0.000001 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 58 0.0082 0.000141 0.0131 0.000226 0.07 0.11 |
| set_parent_processor_ids() 32 0.0015 0.000047 0.0015 0.000047 0.01 0.01 |
| |
| PatchRecoveryErrorEstimator |
| estimate_error() 120 1.5196 0.012663 6.1213 0.051011 12.42 50.04 |
| |
| PetscLinearSolver |
| solve() 14 0.3192 0.022799 0.3192 0.022799 2.61 2.61 |
| |
| ProjectVector |
| operator() 10 0.0115 0.001148 0.1799 0.017989 0.09 1.47 |
| |
| System |
| project_vector() 10 0.0081 0.000807 0.2667 0.026675 0.07 2.18 |
| |
| WeightedPatchRecoveryErrorEstimator |
| estimate_error() 40 0.7117 0.017793 3.1254 0.078134 5.82 25.55 |
| |
| XdrIO |
| read() 1 0.0015 0.001522 0.0016 0.001599 0.01 0.01 |
---------------------------------------------------------------------------------------------------------------------
| Totals: 250642 12.2328 100.00 |
---------------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adjoints_ex3:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************