Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

fix for initing / seting up preconditioners with adaptivity #38

Merged
merged 25 commits into from

5 participants

@friedmud
Owner

This changes the way Preconditioners work. If people have custom Preconditioners they may need to update their objects to implement setup()

If you want to wait until after the 0.9 release to apply this, that's ok with me.

jwpeterson added some commits
@jwpeterson jwpeterson Changes in adaptivity_ex2.C for -Wshadow.
Issue #24.
b4da774
@jwpeterson jwpeterson Changes in adaptivity_ex3.C for -Wshadow.
Issue #24.
e17c79e
@jwpeterson jwpeterson Changes in adaptivity_ex5.C for -Wshadow.
Issue #24.
729541a
@jwpeterson jwpeterson Changes in adjoints_ex1 for -Wshadow.
Issue #24.
b657617
@jwpeterson jwpeterson Changes in adjoints_ex2 for -Wshadow.
Issue #24.
e6f444c
@jwpeterson jwpeterson Changes in adjoints_ex3 for -Wshadow.
Issue #24.
cccdb62
@jwpeterson jwpeterson Changes in adjoints_ex4 for -Wshadow.
Issue #24.
2b9c717
@jwpeterson jwpeterson Don't include deprecated o_string_stream header. ba6afdd
@jwpeterson jwpeterson Fixes for unused variable warnings. 0dc8e2e
@jwpeterson jwpeterson Changes in adjoints_ex5 for -Wshadow.
Issue #24.
f957dea
@jwpeterson jwpeterson Changes in fem_system_ex1 for -Wshadow.
Issue #24.
6ccd84e
@jwpeterson jwpeterson Changes in fem_system_ex2 for -Wshadow.
Issue #24.
85a7e97
@jwpeterson jwpeterson Changes in miscellaneous_ex3 for -Wshadow.
Issue #24.
eb616e8
@jwpeterson jwpeterson Changes in type_vector.h for -Wshadow.
Issue #24.
fe3039d
@jwpeterson jwpeterson Changes in miscellaneous_ex4.C for -Wshadow.
Issue #24.
1b4c2a7
@jwpeterson jwpeterson Changes in miscellaneous_ex7.C for -Wshadow.
Issue #24.
2375c2c
@jwpeterson jwpeterson Changes in reduced_basis_ex1 for -Wshadow.
Issue #24.
0392b90
@jwpeterson jwpeterson Changes in reduced_basis_ex3 for -Wshadow.
Issue #24.
95d8003
@jwpeterson jwpeterson Changes in reduced_basis_ex4 for -Wshadow.
Issue #24.
1a5db52
@jwpeterson jwpeterson Changes in reduced_basis_ex5 for -Wshadow.
Issue #24.
3847f57
@jwpeterson jwpeterson Changes in transient_ex1 for -Wshadow.
Issue #24.
5bcca40
@jwpeterson jwpeterson Changes in vector_fe_ex2 for -Wshadow.
Issue #24.
5e0f961
@jwpeterson jwpeterson Changes in vector_fe_ex3 for -Wshadow.
Issue #24.
14f5c71
@jwpeterson jwpeterson Changes in reduced_basis_ex6 for -Wshadow.
Issue #24.
d99e3fc
@benkirk
Owner

I'm probably being paranoid, but would you mind changing the " out << " to " err << " ?

Otherwise I agree with @roystgnr, this is worth getting in 0.9.0.

@friedmud
Owner

Ok I'll fix that (and I saw that I wasn't checking the return of PCDestroy() so I'll fix that up as well) and then I'll push it in a bit later tonight.

@friedmud
Owner

Ok I pushed this manually... The hash is 213e204

BTW - if we fetch and rebase on master and push manually using

git push upstream HEAD:master

We won't have merge commits on Master. Just another friendly reminder......

@friedmud friedmud closed this
@roystgnr
Owner
@benkirk
Owner
@friedmud
Owner

;-)

I still really hate merge commits.

At best they're confusing (they pull old changes forward... making it look like a change is happening again)

At worst they screw with bisections... and with backtracking and undoing changes in general.

@roystgnr
Owner
@benkirk
Owner

I'm seeing this

  CXX    src/numerics/libmesh_opt_la-petsc_vector.lo
src/numerics/petsc_preconditioner.C: In member function 'virtual void libMesh::PetscPreconditioner<T>::init()':
src/numerics/petsc_preconditioner.C:74:32: error: cannot convert '_p_PC**' to 'PC {aka _p_PC*}' for argument '1' to 'PetscErrorCode PCDestroy(PC)'
make[1]: *** [src/numerics/libmesh_opt_la-petsc_preconditioner.lo] Error 1

on ubuntu-LTS with PETSc-3.1.

I loop over petsc 3.0, 3.1, 3.2, & 3.3 in a different test setup, but that will run tonight.

I'll see what else pops up.

@benkirk benkirk reopened this
@benkirk benkirk merged commit 213e204 into libMesh:master
@friedmud
Owner

Ben - I don't understand what you did here. Why did you open this and remerge it? I'm thoroughly confused.

As for the petsc issues... it might be an issue with 3.1 vs 3.3. I'll look into it in the morning

@benkirk
Owner
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Feb 5, 2013
  1. @jwpeterson
  2. @jwpeterson
  3. @jwpeterson
  4. @jwpeterson
  5. @jwpeterson
  6. @jwpeterson
  7. @jwpeterson
  8. @jwpeterson
  9. @jwpeterson
  10. @jwpeterson
  11. @jwpeterson
  12. @jwpeterson
  13. @jwpeterson
  14. @jwpeterson
  15. @jwpeterson
  16. @jwpeterson
  17. @jwpeterson
  18. @jwpeterson
  19. @jwpeterson
  20. @jwpeterson
  21. @jwpeterson
  22. @jwpeterson
  23. @jwpeterson
  24. @jwpeterson
  25. @friedmud
This page is out of date. Refresh to see the latest.
Showing with 166 additions and 138 deletions.
  1. +0 −3  examples/adaptivity/adaptivity_ex2/adaptivity_ex2.C
  2. +6 −6 examples/adaptivity/adaptivity_ex3/adaptivity_ex3.C
  3. +2 −9 examples/adaptivity/adaptivity_ex5/adaptivity_ex5.C
  4. +3 −3 examples/adjoints/adjoints_ex1/L-shaped.h
  5. +3 −3 examples/adjoints/adjoints_ex2/L-shaped.h
  6. +6 −7 examples/adjoints/adjoints_ex2/adjoints_ex2.C
  7. +3 −3 examples/adjoints/adjoints_ex3/coupled_system.h
  8. +3 −3 examples/adjoints/adjoints_ex4/L-shaped.h
  9. +1 −4 examples/adjoints/adjoints_ex5/adjoints_ex5.C
  10. +4 −4 examples/adjoints/adjoints_ex5/factoryfunction.C
  11. +14 −14 examples/adjoints/adjoints_ex5/femparameters.C
  12. +6 −6 examples/adjoints/adjoints_ex5/heatsystem.C
  13. +3 −3 examples/adjoints/adjoints_ex5/heatsystem.h
  14. +3 −3 examples/fem_system/fem_system_ex1/naviersystem.h
  15. +2 −2 examples/fem_system/fem_system_ex2/nonlinear_neohooke_cc.h
  16. +3 −3 examples/fem_system/fem_system_ex2/solid_system.C
  17. +2 −2 examples/miscellaneous/miscellaneous_ex3/miscellaneous_ex3.C
  18. +0 −1  examples/miscellaneous/miscellaneous_ex4/miscellaneous_ex4.C
  19. +3 −3 examples/miscellaneous/miscellaneous_ex7/biharmonic.C
  20. +3 −3 examples/miscellaneous/miscellaneous_ex7/biharmonic_jr.C
  21. +3 −3 examples/reduced_basis/reduced_basis_ex1/rb_classes.h
  22. +3 −3 examples/reduced_basis/reduced_basis_ex3/rb_classes.h
  23. +3 −3 examples/reduced_basis/reduced_basis_ex4/eim_classes.h
  24. +3 −3 examples/reduced_basis/reduced_basis_ex4/rb_classes.h
  25. +3 −3 examples/reduced_basis/reduced_basis_ex5/rb_classes.h
  26. +0 −2  examples/reduced_basis/reduced_basis_ex5/reduced_basis_ex5.C
  27. +4 −4 examples/reduced_basis/reduced_basis_ex6/eim_classes.h
  28. +3 −3 examples/reduced_basis/reduced_basis_ex6/rb_classes.h
  29. +0 −3  examples/transient/transient_ex1/transient_ex1.C
  30. +3 −3 examples/vector_fe/vector_fe_ex2/laplace_system.C
  31. +4 −4 examples/vector_fe/vector_fe_ex2/solution_function.h
  32. +3 −3 examples/vector_fe/vector_fe_ex3/curl_curl_system.C
  33. +4 −4 examples/vector_fe/vector_fe_ex3/solution_function.h
  34. +9 −0 include/numerics/preconditioner.h
  35. +4 −4 include/numerics/type_vector.h
  36. +6 −3 src/numerics/petsc_preconditioner.C
  37. +34 −3 src/solvers/petsc_linear_solver.C
  38. +7 −2 src/solvers/petsc_nonlinear_solver.C
View
3  examples/adaptivity/adaptivity_ex2/adaptivity_ex2.C
@@ -440,9 +440,6 @@ int main (int argc, char** argv)
if(!read_solution)
{
// Print out the H1 norm of the saved solution, for verification purposes:
- TransientLinearImplicitSystem& system =
- equation_systems.get_system<TransientLinearImplicitSystem>
- ("Convection-Diffusion");
Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "Final H1 norm = " << H1norm << std::endl << std::endl;
View
12 examples/adaptivity/adaptivity_ex3/adaptivity_ex3.C
@@ -578,7 +578,7 @@ void assemble_laplace(EquationSystems& es,
const MeshBase& mesh = es.get_mesh();
// The dimension that we are running
- const unsigned int dim = mesh.mesh_dimension();
+ const unsigned int mesh_dim = mesh.mesh_dimension();
// Get a reference to the LinearImplicitSystem we are solving
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Laplace");
@@ -597,12 +597,12 @@ void assemble_laplace(EquationSystems& es,
// \p FEBase::build() member dynamically creates memory we will
// store the object as an \p AutoPtr<FEBase>. This can be thought
// of as a pointer that will clean up after itself.
- AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
- AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
+ AutoPtr<FEBase> fe (FEBase::build(mesh_dim, fe_type));
+ AutoPtr<FEBase> fe_face (FEBase::build(mesh_dim, fe_type));
// Quadrature rules for numerical integration.
- AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(dim));
- AutoPtr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
+ AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(mesh_dim));
+ AutoPtr<QBase> qface(fe_type.default_quadrature_rule(mesh_dim-1));
// Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (qrule.get());
@@ -709,7 +709,7 @@ void assemble_laplace(EquationSystems& es,
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
// We need a forcing function to make the 1D case interesting
- if (dim == 1)
+ if (mesh_dim == 1)
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
{
Real x = q_point[qp](0);
View
11 examples/adaptivity/adaptivity_ex5/adaptivity_ex5.C
@@ -415,9 +415,6 @@ int main (int argc, char** argv)
// vector assignment. Since only \p TransientLinearImplicitSystems
// (and systems derived from them) contain old solutions
// we need to specify the system type when we ask for it.
- TransientLinearImplicitSystem & system =
- equation_systems.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
-
*system.old_local_solution = *system.current_local_solution;
// The number of refinement steps per time step.
@@ -432,8 +429,7 @@ int main (int argc, char** argv)
system.solve();
// Print out the H1 norm, for verification purposes:
- Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
-
+ H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "H1 norm = " << H1norm << std::endl;
// Possibly refine the mesh
@@ -515,10 +511,7 @@ int main (int argc, char** argv)
if(!read_solution)
{
// Print out the H1 norm of the saved solution, for verification purposes:
- TransientLinearImplicitSystem& system =
- equation_systems.get_system<TransientLinearImplicitSystem>
- ("Convection-Diffusion");
- Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
+ H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "Final H1 norm = " << H1norm << std::endl << std::endl;
View
6 examples/adjoints/adjoints_ex1/L-shaped.h
@@ -13,9 +13,9 @@ class LaplaceSystem : public FEMSystem
public:
// Constructor
LaplaceSystem(EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : FEMSystem(es, name_in, number_in),
_fe_family("LAGRANGE"), _fe_order(1),
_analytic_jacobians(true) { qoi.resize(2); }
View
6 examples/adjoints/adjoints_ex2/L-shaped.h
@@ -14,9 +14,9 @@ class LaplaceSystem : public FEMSystem
public:
// Constructor
LaplaceSystem(EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : FEMSystem(es, name_in, number_in),
_fe_family("LAGRANGE"), _fe_order(1),
_analytic_jacobians(true) { }
View
13 examples/adjoints/adjoints_ex2/adjoints_ex2.C
@@ -263,7 +263,6 @@ int main (int argc, char** argv)
}
}
GetPot infile("general.in");
- GetPot infile_l_shaped("l-shaped.in");
// Read in parameters from the input file
FEMParameters param;
@@ -364,12 +363,12 @@ int main (int argc, char** argv)
// 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);
- GetPot infile("l-shaped.in");
+ GetPot infile_l_shaped("l-shaped.in");
Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
- Number sensitivity_QoI_0_0_exact = infile("sensitivity_0_0", 0.0);
+ Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
- Number sensitivity_QoI_0_1_exact = infile("sensitivity_0_1", 0.0);
+ Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
<< " active elements and "
@@ -482,12 +481,12 @@ int main (int argc, char** argv)
// 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);
- GetPot infile("l-shaped.in");
+ GetPot infile_l_shaped("l-shaped.in");
Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
- Number sensitivity_QoI_0_0_exact = infile("sensitivity_0_0", 0.0);
+ Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
- Number sensitivity_QoI_0_1_exact = infile("sensitivity_0_1", 0.0);
+ Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
<< " active elements and "
View
6 examples/adjoints/adjoints_ex3/coupled_system.h
@@ -34,9 +34,9 @@ class CoupledSystem : public FEMSystem
public:
// Constructor
CoupledSystem(EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number), Peclet(1.) {qoi.resize(1);}
+ 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
View
6 examples/adjoints/adjoints_ex4/L-shaped.h
@@ -13,9 +13,9 @@ class LaplaceSystem : public FEMSystem
public:
// Constructor
LaplaceSystem(EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : FEMSystem(es, name_in, number_in),
_fe_family("LAGRANGE"), _fe_order(1),
_analytic_jacobians(true) { qoi.resize(2); }
View
5 examples/adjoints/adjoints_ex5/adjoints_ex5.C
@@ -96,13 +96,10 @@
#include "libmesh/solution_history.h"
#include "libmesh/memory_solution_history.h"
-// Some (older) compilers do not offer full stream
-// functionality, OStringStream works around this.
-#include "libmesh/o_string_stream.h"
-
// C++ includes
#include <iostream>
#include <sys/time.h>
+#include <iomanip>
void write_output(EquationSystems &es,
unsigned int a_step, // The adaptive step count
View
8 examples/adjoints/adjoints_ex5/factoryfunction.C
@@ -28,14 +28,14 @@ using namespace libMesh;
class ExampleOneFunction : public FunctionBase<Number>
{
- virtual Number operator() (const Point& p,
- const Real time = 0)
+ virtual Number operator() (const Point& /*p*/,
+ const Real /*time*/)
{
return 1;
}
- virtual void operator() (const Point& p,
- const Real time,
+ virtual void operator() (const Point& /*p*/,
+ const Real /*time*/,
DenseVector<Number>& output)
{
for (unsigned int i=0; i != output.size(); ++i)
View
28 examples/adjoints/adjoints_ex5/femparameters.C
@@ -282,17 +282,17 @@ void FEMParameters::read(GetPot &input)
{
unsigned int myboundary =
input("periodic_boundaries", -1, i);
- unsigned int pairedboundary = 0;
+ // unsigned int pairedboundary = 0;
RealVectorValue translation_vector;
if (dimension == 2)
switch (myboundary)
{
case 0:
- pairedboundary = 2;
+ // pairedboundary = 2;
translation_vector = RealVectorValue(0., domain_edge_length);
break;
case 1:
- pairedboundary = 3;
+ // pairedboundary = 3;
translation_vector = RealVectorValue(-domain_edge_width, 0);
break;
default:
@@ -304,15 +304,15 @@ void FEMParameters::read(GetPot &input)
switch (myboundary)
{
case 0:
- pairedboundary = 5;
+ // pairedboundary = 5;
translation_vector = RealVectorValue(0., 0., domain_edge_height);
break;
case 1:
- pairedboundary = 3;
+ // pairedboundary = 3;
translation_vector = RealVectorValue(0., domain_edge_length, 0.);
break;
case 2:
- pairedboundary = 4;
+ // pairedboundary = 4;
translation_vector = RealVectorValue(-domain_edge_width, 0., 0.);
break;
default:
@@ -389,11 +389,11 @@ void FEMParameters::read(GetPot &input)
const std::string variable_set =
input("dirichlet_condition_variables", empty_string, i);
- for (unsigned int i=0; i != variable_set.size(); ++i)
+ for (unsigned int j=0; j != variable_set.size(); ++j)
{
- if (variable_set[i] == '1')
- dirichlet_condition_variables[func_boundary].push_back(i);
- else if (variable_set[i] != '0')
+ if (variable_set[j] == '1')
+ dirichlet_condition_variables[func_boundary].push_back(j);
+ else if (variable_set[j] != '0')
{
libMesh::out << "Unable to understand Dirichlet variable set"
<< variable_set << std::endl;
@@ -460,11 +460,11 @@ void FEMParameters::read(GetPot &input)
const std::string variable_set =
input("neumann_condition_variables", empty_string, i);
- for (unsigned int i=0; i != variable_set.size(); ++i)
+ for (unsigned int j=0; j != variable_set.size(); ++j)
{
- if (variable_set[i] == '1')
- neumann_condition_variables[func_boundary].push_back(i);
- else if (variable_set[i] != '0')
+ if (variable_set[j] == '1')
+ neumann_condition_variables[func_boundary].push_back(j);
+ else if (variable_set[j] != '0')
{
libMesh::out << "Unable to understand Neumann variable set"
<< variable_set << std::endl;
View
12 examples/adjoints/adjoints_ex5/heatsystem.C
@@ -157,15 +157,15 @@ bool HeatSystem::element_time_derivative (bool request_jacobian,
}
// Perturb and accumulate dual weighted residuals
-void HeatSystem::perturb_accumulate_residuals(const ParameterVector& parameters)
+void HeatSystem::perturb_accumulate_residuals(const ParameterVector& parameters_in)
{
- const unsigned int Np = parameters.size();
+ const unsigned int Np = parameters_in.size();
for (unsigned int j=0; j != Np; ++j)
{
- Number old_parameter = *parameters[j];
+ Number old_parameter = *parameters_in[j];
- *parameters[j] = old_parameter - dp;
+ *parameters_in[j] = old_parameter - dp;
this->assembly(true, false);
@@ -177,7 +177,7 @@ void HeatSystem::perturb_accumulate_residuals(const ParameterVector& parameters)
// But since we compute the residual already scaled by dt, there is no need for the * dt
R_minus_dp += -R_minus->dot(this->get_adjoint_solution(0));
- *parameters[j] = old_parameter + dp;
+ *parameters_in[j] = old_parameter + dp;
this->assembly(true, false);
@@ -187,7 +187,7 @@ void HeatSystem::perturb_accumulate_residuals(const ParameterVector& parameters)
R_plus_dp += -R_plus->dot(this->get_adjoint_solution(0));
- *parameters[j] = old_parameter;
+ *parameters_in[j] = old_parameter;
}
}
View
6 examples/adjoints/adjoints_ex5/heatsystem.h
@@ -28,9 +28,9 @@ class HeatSystem : public FEMSystem
public:
// Constructor
HeatSystem(EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : FEMSystem(es, name_in, number_in),
_k(1.0),
_fe_family("LAGRANGE"), _fe_order(1),
_analytic_jacobians(true), R_plus_dp(0.0), R_minus_dp(0.0), dp(1.e-6) { qoi.resize(1); }
View
6 examples/fem_system/fem_system_ex1/naviersystem.h
@@ -28,9 +28,9 @@ class NavierSystem : public FEMSystem
public:
// Constructor
NavierSystem(EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number), Reynolds(1.), application(0) {}
+ const std::string& name_in,
+ const unsigned int number_in)
+ : FEMSystem(es, name_in, number_in), Reynolds(1.), application(0) {}
// System initialization
virtual void init_data ();
View
4 examples/fem_system/fem_system_ex2/nonlinear_neohooke_cc.h
@@ -37,8 +37,8 @@ using namespace libMesh;
class NonlinearNeoHookeCurrentConfig {
public:
NonlinearNeoHookeCurrentConfig(
- const std::vector<std::vector<RealGradient> >& dphi, GetPot& args) :
- dphi(dphi) {
+ const std::vector<std::vector<RealGradient> >& dphi_in, GetPot& args) :
+ dphi(dphi_in) {
E = args("material/neohooke/e_modulus", 10000.0);
nu = args("material/neohooke/nu", 0.3);
}
View
6 examples/fem_system/fem_system_ex2/solid_system.C
@@ -43,9 +43,9 @@
// Bring in everything from the libMesh namespace
using namespace libMesh;
-SolidSystem::SolidSystem(EquationSystems& es, const std::string& name,
- const unsigned int number) :
- FEMSystem(es, name, number) {
+SolidSystem::SolidSystem(EquationSystems& es, const std::string& name_in,
+ const unsigned int number_in) :
+ FEMSystem(es, name_in, number_in) {
// Add a time solver. We are just looking at a steady state problem here.
this->time_solver = AutoPtr<TimeSolver>(new SteadySolver(*this));
View
4 examples/miscellaneous/miscellaneous_ex3/miscellaneous_ex3.C
@@ -83,7 +83,7 @@ const Real sigma = 0.2;
// This function computes the Jacobian K(x)
void compute_jacobian (const NumericVector<Number>& soln,
SparseMatrix<Number>& jacobian,
- NonlinearImplicitSystem& sys)
+ NonlinearImplicitSystem& /*sys*/)
{
// Get a reference to the equation system.
EquationSystems &es = *_equation_system;
@@ -213,7 +213,7 @@ void compute_jacobian (const NumericVector<Number>& soln,
// x is passed in the soln vector
void compute_residual (const NumericVector<Number>& soln,
NumericVector<Number>& residual,
- NonlinearImplicitSystem& sys)
+ NonlinearImplicitSystem& /*sys*/)
{
EquationSystems &es = *_equation_system;
View
1  examples/miscellaneous/miscellaneous_ex4/miscellaneous_ex4.C
@@ -110,7 +110,6 @@ int main (int argc, char** argv)
// Create an equation systems object.
EquationSystems equation_systems (mesh);
- MeshRefinement mesh_refinement (mesh);
MeshTools::Generation::build_square (mesh,
16,
View
6 examples/miscellaneous/miscellaneous_ex7/biharmonic.C
@@ -96,16 +96,16 @@ void Biharmonic::init()
-void Biharmonic::step(const Real& dt)
+void Biharmonic::step(const Real& dt_in)
{
// We need to update the old solution vector.
// The old solution vector will be the current solution vector from the
// previous time step. We use vector assignment. Only \p TransientSystems
// (and systems derived from them) contain old solutions.
- if (dt < 0)
+ if (dt_in < 0)
_dt = _dt0;
else
- _dt = dt;
+ _dt = dt_in;
*(_jr->old_local_solution) = *(_jr->current_local_solution);
View
6 examples/miscellaneous/miscellaneous_ex7/biharmonic_jr.C
@@ -16,9 +16,9 @@
using namespace libMesh;
Biharmonic::JR::JR(EquationSystems& eqSys,
- const std::string& name,
- const unsigned int number) :
- TransientNonlinearImplicitSystem(eqSys,name,number),
+ const std::string& name_in,
+ const unsigned int number_in) :
+ TransientNonlinearImplicitSystem(eqSys,name_in,number_in),
_biharmonic(dynamic_cast<Biharmonic&>(eqSys))
{
// Check that we can actually compute second derivatives
View
6 examples/reduced_basis/reduced_basis_ex1/rb_classes.h
@@ -73,9 +73,9 @@ class SimpleRBConstruction : public RBConstruction
public:
SimpleRBConstruction (EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : Parent(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : Parent(es, name_in, number_in),
dirichlet_bc(AutoPtr<DirichletBoundary>(NULL))
{}
View
6 examples/reduced_basis/reduced_basis_ex3/rb_classes.h
@@ -73,9 +73,9 @@ class SimpleRBConstruction : public TransientRBConstruction
public:
SimpleRBConstruction (EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : Parent(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : Parent(es, name_in, number_in),
dirichlet_bc(AutoPtr<DirichletBoundary>(NULL))
{}
View
6 examples/reduced_basis/reduced_basis_ex4/eim_classes.h
@@ -39,9 +39,9 @@ class SimpleEIMConstruction : public RBEIMConstruction
* Constructor.
*/
SimpleEIMConstruction (EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : Parent(es, name, number)
+ const std::string& name_in,
+ const unsigned int number_in)
+ : Parent(es, name_in, number_in)
{
}
View
6 examples/reduced_basis/reduced_basis_ex4/rb_classes.h
@@ -55,9 +55,9 @@ class SimpleRBConstruction : public RBConstruction
public:
SimpleRBConstruction (EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : Parent(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : Parent(es, name_in, number_in),
dirichlet_bc(AutoPtr<DirichletBoundary>(NULL))
{}
View
6 examples/reduced_basis/reduced_basis_ex5/rb_classes.h
@@ -45,9 +45,9 @@ class ElasticityRBConstruction : public RBConstruction
public:
ElasticityRBConstruction (EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : Parent(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : Parent(es, name_in, number_in),
elasticity_assembly_expansion(*this),
ip_assembly(*this)
{}
View
2  examples/reduced_basis/reduced_basis_ex5/reduced_basis_ex5.C
@@ -208,8 +208,6 @@ int main(int argc, char** argv) {
void scale_mesh_and_plot(EquationSystems& es, const RBParameters& mu, const std::string& filename)
{
- const Real x_scaling = mu.get_value("x_scaling");
-
// Loop over the mesh nodes and move them!
MeshBase& mesh = es.get_mesh();
View
8 examples/reduced_basis/reduced_basis_ex6/eim_classes.h
@@ -46,9 +46,9 @@ class SimpleEIMConstruction : public RBEIMConstruction
* Constructor.
*/
SimpleEIMConstruction (EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : Parent(es, name, number)
+ const std::string& name_in,
+ const unsigned int number_in)
+ : Parent(es, name_in, number_in)
{
}
@@ -93,4 +93,4 @@ class SimpleEIMConstruction : public RBEIMConstruction
};
-#endif
+#endif
View
6 examples/reduced_basis/reduced_basis_ex6/rb_classes.h
@@ -74,9 +74,9 @@ class SimpleRBConstruction : public RBConstruction
public:
SimpleRBConstruction (EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : Parent(es, name, number),
+ const std::string& name_in,
+ const unsigned int number_in)
+ : Parent(es, name_in, number_in),
ex6_assembly_expansion(*this),
dirichlet_bc(AutoPtr<DirichletBoundary>(NULL))
{}
View
3  examples/transient/transient_ex1/transient_ex1.C
@@ -218,9 +218,6 @@ int main (int argc, char** argv)
// vector assignment. Since only \p TransientSystems
// (and systems derived from them) contain old solutions
// we need to specify the system type when we ask for it.
- TransientLinearImplicitSystem& system =
- equation_systems.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
-
*system.old_local_solution = *system.current_local_solution;
// Assemble & solve the linear system
View
6 examples/vector_fe/vector_fe_ex2/laplace_system.C
@@ -33,9 +33,9 @@
using namespace libMesh;
LaplaceSystem::LaplaceSystem( EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number)
+ const std::string& name_in,
+ const unsigned int number_in)
+ : FEMSystem(es, name_in, number_in)
{
return;
}
View
8 examples/vector_fe/vector_fe_ex2/solution_function.h
@@ -47,11 +47,11 @@ class SolutionFunction : public FunctionBase<Number>
output(_u_var+2) = soln( 2, x, y, z );
}
- virtual Number component( unsigned int component, const Point& p,
+ virtual Number component( unsigned int component_in, const Point& p,
const Real )
{
const Real x=p(0), y=p(1), z=p(2);
- return soln( component, x, y, z );
+ return soln( component_in, x, y, z );
}
virtual AutoPtr<FunctionBase<Number> > clone() const
@@ -87,11 +87,11 @@ class SolutionGradient : public FunctionBase<Gradient>
output(_u_var+2) = soln( 2, x, y, z );
}
- virtual Gradient component( unsigned int component, const Point& p,
+ virtual Gradient component( unsigned int component_in, const Point& p,
const Real )
{
const Real x=p(0), y=p(1), z=p(2);
- return soln( component, x, y, z );
+ return soln( component_in, x, y, z );
}
virtual AutoPtr<FunctionBase<Gradient> > clone() const
View
6 examples/vector_fe/vector_fe_ex3/curl_curl_system.C
@@ -33,9 +33,9 @@
using namespace libMesh;
CurlCurlSystem::CurlCurlSystem( EquationSystems& es,
- const std::string& name,
- const unsigned int number)
- : FEMSystem(es, name, number)
+ const std::string& name_in,
+ const unsigned int number_in)
+ : FEMSystem(es, name_in, number_in)
{
return;
}
View
8 examples/vector_fe/vector_fe_ex3/solution_function.h
@@ -46,11 +46,11 @@ class SolutionFunction : public FunctionBase<Number>
output(_u_var+1) = soln( x, y )(1);
}
- virtual Number component( unsigned int component, const Point& p,
+ virtual Number component( unsigned int component_in, const Point& p,
const Real )
{
const Real x=p(0), y=p(1);
- return soln( x, y )(component);
+ return soln( x, y )(component_in);
}
virtual AutoPtr<FunctionBase<Number> > clone() const
@@ -82,11 +82,11 @@ class SolutionGradient : public FunctionBase<Gradient>
output(_u_var+1) = soln.grad( x, y ).row(1);
}
- virtual Gradient component( unsigned int component, const Point& p,
+ virtual Gradient component( unsigned int component_in, const Point& p,
const Real )
{
const Real x=p(0), y=p(1);
- return soln.grad( x, y ).row(component);
+ return soln.grad( x, y ).row(component_in);
}
virtual AutoPtr<FunctionBase<Gradient> > clone() const
View
9 include/numerics/preconditioner.h
@@ -98,10 +98,19 @@ class Preconditioner : public ReferenceCountedObject<Preconditioner<T> >
/**
* Initialize data structures if not done so already.
+ *
+ * This MUST be called before the preconditioning object is used.
*/
virtual void init () {}
/**
+ * This is called every time the "operator might have changed".
+ *
+ * This is essentially where you need to fill in your preconditioning matrix.
+ */
+ virtual void setup () {}
+
+ /**
* Sets the matrix P to be preconditioned.
*/
void set_matrix(SparseMatrix<Number> & mat);
View
8 include/numerics/type_vector.h
@@ -651,19 +651,19 @@ typename boostcopy::enable_if_c<
TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
TypeVector<T>::operator * (const Scalar factor) const
{
- typedef typename CompareTypes<T, Scalar>::supertype TS;
+ typedef typename CompareTypes<T, Scalar>::supertype SuperType;
#if LIBMESH_DIM == 1
- return TypeVector<TS>(_coords[0]*factor);
+ return TypeVector<SuperType>(_coords[0]*factor);
#endif
#if LIBMESH_DIM == 2
- return TypeVector<TS>(_coords[0]*factor,
+ return TypeVector<SuperType>(_coords[0]*factor,
_coords[1]*factor);
#endif
#if LIBMESH_DIM == 3
- return TypeVector<TS>(_coords[0]*factor,
+ return TypeVector<SuperType>(_coords[0]*factor,
_coords[1]*factor,
_coords[2]*factor);
#endif
View
9 src/numerics/petsc_preconditioner.C
@@ -68,13 +68,16 @@ void PetscPreconditioner<T>::init ()
// Clear the preconditioner in case it has been created in the past
if (!this->_is_initialized)
{
- // Create the preconditioning object
- if (!_pc)
+ // Should probably use PCReset(), but it's not working at the moment so we'll destroy instead
+ if (_pc)
{
- int ierr = PCCreate(libMesh::COMM_WORLD,&_pc);
+ int ierr = PCDestroy(&_pc);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
+ int ierr = PCCreate(libMesh::COMM_WORLD,&_pc);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
PetscMatrix<T> * pmatrix = libmesh_cast_ptr<PetscMatrix<T>*, SparseMatrix<T> >(this->_matrix);
_mat = pmatrix->mat();
View
37 src/solvers/petsc_linear_solver.C
@@ -46,7 +46,14 @@ extern "C"
PetscErrorCode __libmesh_petsc_preconditioner_setup (void * ctx)
{
Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
- preconditioner->init();
+
+ if(!preconditioner->initialized())
+ {
+ err<<"Preconditioner not initialized! Make sure you call init() before solve!"<<std::endl;
+ libmesh_error();
+ }
+
+ preconditioner->setup();
return 0;
}
@@ -69,7 +76,14 @@ extern "C"
void *ctx;
PetscErrorCode ierr = PCShellGetContext(pc,&ctx);CHKERRQ(ierr);
Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
- preconditioner->init();
+
+ if(!preconditioner->initialized())
+ {
+ err<<"Preconditioner not initialized! Make sure you call init() before solve!"<<std::endl;
+ libmesh_error();
+ }
+
+ preconditioner->setup();
return 0;
}
@@ -249,6 +263,7 @@ void PetscLinearSolver<T>::init ()
//If there is a preconditioner object we need to set the internal setup and apply routines
if(this->_preconditioner)
{
+ this->_preconditioner->init();
PCShellSetContext(_pc,(void*)this->_preconditioner);
PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
@@ -366,6 +381,7 @@ void PetscLinearSolver<T>::init ( PetscMatrix<T>* matrix )
if(this->_preconditioner)
{
this->_preconditioner->set_matrix(*matrix);
+ this->_preconditioner->init();
PCShellSetContext(_pc,(void*)this->_preconditioner);
PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
@@ -648,6 +664,7 @@ PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in,
{
subprecond_matrix = new PetscMatrix<Number>(subprecond);
this->_preconditioner->set_matrix(*subprecond_matrix);
+ this->_preconditioner->init();
}
}
else
@@ -657,7 +674,10 @@ PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in,
CHKERRABORT(libMesh::COMM_WORLD,ierr);
if(this->_preconditioner)
+ {
this->_preconditioner->set_matrix(matrix_in);
+ this->_preconditioner->init();
+ }
}
// Set the tolerances for the iterative solver. Use the user-supplied
@@ -718,6 +738,7 @@ PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in,
/* Before we delete subprecond_matrix, we should give the
_preconditioner a different matrix. */
this->_preconditioner->set_matrix(matrix_in);
+ this->_preconditioner->init();
delete subprecond_matrix;
subprecond_matrix = NULL;
}
@@ -980,6 +1001,7 @@ PetscLinearSolver<T>::adjoint_solve (SparseMatrix<T>& matrix_in,
{
subprecond_matrix = new PetscMatrix<Number>(subprecond);
this->_preconditioner->set_matrix(*subprecond_matrix);
+ this->_preconditioner->init();
}
}
else
@@ -989,7 +1011,10 @@ PetscLinearSolver<T>::adjoint_solve (SparseMatrix<T>& matrix_in,
CHKERRABORT(libMesh::COMM_WORLD,ierr);
if(this->_preconditioner)
+ {
this->_preconditioner->set_matrix(matrix_in);
+ this->_preconditioner->init();
+ }
}
// Set the tolerances for the iterative solver. Use the user-supplied
@@ -1050,6 +1075,7 @@ PetscLinearSolver<T>::adjoint_solve (SparseMatrix<T>& matrix_in,
/* Before we delete subprecond_matrix, we should give the
_preconditioner a different matrix. */
this->_preconditioner->set_matrix(matrix_in);
+ this->_preconditioner->init();
delete subprecond_matrix;
subprecond_matrix = NULL;
}
@@ -1552,6 +1578,7 @@ PetscLinearSolver<T>::solve (const ShellMatrix<T>& shell_matrix,
{
subprecond_matrix = new PetscMatrix<Number>(subprecond);
this->_preconditioner->set_matrix(*subprecond_matrix);
+ this->_preconditioner->init();
}
}
else
@@ -1561,7 +1588,10 @@ PetscLinearSolver<T>::solve (const ShellMatrix<T>& shell_matrix,
CHKERRABORT(libMesh::COMM_WORLD,ierr);
if(this->_preconditioner)
+ {
this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number>&>(precond_matrix));
+ this->_preconditioner->init();
+ }
}
// Set the tolerances for the iterative solver. Use the user-supplied
@@ -1622,6 +1652,7 @@ PetscLinearSolver<T>::solve (const ShellMatrix<T>& shell_matrix,
/* Before we delete subprecond_matrix, we should give the
_preconditioner a different matrix. */
this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number>&>(precond_matrix));
+ this->_preconditioner->init();
delete subprecond_matrix;
subprecond_matrix = NULL;
}
@@ -1760,7 +1791,7 @@ void PetscLinearSolver<T>::set_petsc_solver_type()
#else
ierr = KSPSetType (_ksp, (char*) KSPCHEBYSHEV); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
#endif
-
+
default:
libMesh::err << "ERROR: Unsupported PETSC Solver: "
View
9 src/solvers/petsc_nonlinear_solver.C
@@ -329,6 +329,8 @@ void PetscNonlinearSolver<T>::init ()
ierr = KSPGetPC(ksp,&pc);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ this->_preconditioner->init();
+
PCSetType(pc, PCSHELL);
PCShellSetContext(pc,(void*)this->_preconditioner);
@@ -462,11 +464,11 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T>& jac_in, // System Jacobian Ma
{
MatNullSpace msp = PETSC_NULL;
this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, &msp);
-
+
if(msp) {
ierr = MatSetNearNullSpace(jac->mat(), msp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
+
ierr = MatNullSpaceDestroy(&msp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
@@ -497,7 +499,10 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T>& jac_in, // System Jacobian Ma
//Set the preconditioning matrix
if(this->_preconditioner)
+ {
this->_preconditioner->set_matrix(jac_in);
+ this->_preconditioner->init();
+ }
// ierr = KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
// CHKERRABORT(libMesh::COMM_WORLD,ierr);
Something went wrong with that request. Please try again.