bgranzow / goal Goto Github PK

View Code? Open in Web Editor NEW

5.0 5.0 2.0 685 KB

Goal Oriented Adaptive Lagrangian Mechanics

License: Other

CMake 9.65% Shell 0.47% C++ 88.61% Makefile 1.27%

adaptation adjoint c-plus-plus error-estimation finite-elements

goal's People

Contributors

Stargazers

Watchers

Forkers

mdipalma 2012zgzyy

goal's Issues

Make a logo

For fun.

Don't place a one on the diagonal for Dirichlet BCs

Rather than zeroing the row corresponding to a Dirichlet BC and placing a one on the diagonal, the diagonal entry could be left as is, and the right hand side scaled by the Jacobian matrix diagonal value. This will likely lead to better conditioned systems.

P-norm functional QoI

elasticity: implement p-norm stress QoI

Implement a QoI of the form:

J(u) = ( \int_{\Omega} | \sigma |_F^p ) ^(1/p)

Get to the bottom of pesky MueLu output

Using the MueLu method: CreateTpetraPreconditioner
with the optional argument inCoords, and specifying the appropriate number of PDE
equations in the input parameterlist inParamList with the parameter number of equations yields the following warning:

******* WARNING *******
Setting matrix block size to 2 (value of the parameter in the list) instead of 1 (provided matrix).
You may want to check "number of equations" (or "PDE equations" for factory style list) parameter.

Omitting the number of equations throws an assertion, presumably because the coordinate vector is assumed to have a 1-1 mapping with the DOF map.

Throw test that evaluated to true: (static_cast<size_t> (X_j_av.size ()) < lclNumRows)

#0  __cxa_throw () at ../../../../gcc-4.9.2/libstdc++-v3/libsupc++/eh_throw.cc:62
#1  0x0000000001c682ba in Tpetra::MultiVector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false>::MultiVector(Teuchos::RCP<Tpetra::Map<int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > const> const&, Teuchos::ArrayView<Teuchos::ArrayView<double const> const> const&, unsigned long) ()
#2  0x000000000120f025 in Xpetra::MultiVectorFactory<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >::Build(Teuchos::RCP<Xpetra::Map<int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > const> const&, Teuchos::ArrayView<Teuchos::ArrayView<double const> const> const&, unsigned long) ()
#3  0x000000000122fe63 in MueLu::Hierarchy<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >::ReplaceCoordinateMap(MueLu::Level&) ()
#4  0x0000000001231821 in MueLu::Hierarchy<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >::Setup(int, Teuchos::RCP<MueLu::FactoryManagerBase const>, Teuchos::RCP<MueLu::FactoryManagerBase const>, Teuchos::RCP<MueLu::FactoryManagerBase const>) ()
#5  0x0000000000ee28ef in MueLu::HierarchyManager<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >::SetupHierarchy(MueLu::Hierarchy<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >&) const ()
#6  0x0000000000ec0895 in Teuchos::RCP<MueLu::Hierarchy<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > > MueLu::CreateXpetraPreconditioner<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >(Teuchos::RCP<Xpetra::Matrix<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > >, Teuchos::ParameterList const&, Teuchos::RCP<Xpetra::MultiVector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > >, Teuchos::RCP<Xpetra::MultiVector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > >) ()
#7  0x0000000000ec1a69 in Teuchos::RCP<MueLu::TpetraOperator<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > > MueLu::CreateTpetraPreconditioner<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> >(Teuchos::RCP<Tpetra::Operator<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace> > > const&, Teuchos::ParameterList&, Teuchos::RCP<Tpetra::MultiVector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>::classic> > const&, Teuchos::RCP<Tpetra::MultiVector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>::classic> > const&) ()
#8  0x0000000000e0a889 in goal::build_multigrid_solver(Teuchos::RCP<Teuchos::ParameterList const>, Teuchos::RCP<goal::Indexer>, Teuchos::RCP<Tpetra::CrsMatrix<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<Tpetra::Vector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<Tpetra::Vector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, int) () at /lore/granzb/goal/src/goal_linear_solvers.cpp:109
#9  0x0000000000e0b407 in goal::solve_multigrid_gmres(Teuchos::RCP<Teuchos::ParameterList const>, Teuchos::RCP<Tpetra::CrsMatrix<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<Tpetra::Vector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<Tpetra::Vector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<goal::Indexer>) () at /lore/granzb/goal/src/goal_linear_solvers.cpp:191
#10 0x0000000000e0bba7 in goal::solve_linear_system(Teuchos::RCP<Teuchos::ParameterList const>, Teuchos::RCP<Tpetra::CrsMatrix<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<Tpetra::Vector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<Tpetra::Vector<double, int, long long, Kokkos::Compat::KokkosDeviceWrapperNode<Kokkos::Serial, Kokkos::HostSpace>, false> >, Teuchos::RCP<goal::Indexer>) () at /lore/granzb/goal/src/goal_linear_solvers.cpp:213
#11 0x0000000000d9a85c in elast::Solver::solve_primal() () at /lore/granzb/goal/miniapps/elasticity/elast_solver.cpp:98
#12 0x0000000000d9b149 in elast::Solver::solve_primal_only() () at /lore/granzb/goal/miniapps/elasticity/elast_solver.cpp:135
#13 0x0000000000d9f2b1 in main () at /lore/granzb/goal/miniapps/elasticity/elast_solver.cpp:164

Why is this happening? I tried to follow very closely (read exactly) what was done in Albany's CTM, where I did not see this issue.

@ibaned, you mentioned you had seen this in the wild before? Any ideas?

elasticity: von mises stress evaluator

Implement iso size field

Based on a single scalar error field.

poisson: implement error estimation

Localize error contributions as:
\mathcal{E}_i = | ( \nabla ( (z^h - z^H) \psi_i), u^H) - ( (z^h-z^H) \psi_i, u^H) |

Update BC documentation for zeroing rows

Now that a 1 is no longer placed on the diagonal. (rather the diagonal entry is left as is).

Documentation page explaining graph-based assembly

A rundown of how Phalanx works with a simple graph example.

Make Indexer an abstract base class

Issue [#24] is making it clear that it would be convenient to have classes like:
StridedIndexer, BlockIndexer, and MixedIndexer all inherit from an abstract Indexer base class. This will probably be painful. Hopefully it's worth it.

Documentation page for mesh generation

Describe the process from:
gmodel -> gmsh -> PUMI

Implement condensation of columns for Dirichlet BCs

For every current envisioned application, the Tpetra CrsGraph (describing the Jacobian matrix sparsity pattern) will be symmetric. This gives a convenient way to zero columns corresponding to prescribed Dirichlet BCs (condensation) which will likely lead to improved linear solver performance. A good starting place to look is here:

https://github.com/dealii/dealii/blob/master/source/numerics/matrix_tools.cc#L168

Fix memory leak in traction bcs

Propogate Goal project version number into Doxyfile from Cmake directly

Implement point-wise scalar quantity of interest

elasticity: implement cauchy stress evaluator

Evaluate the Cauchy with a standard small strain isotropic constitutive relation.

Why do CG and GMRES take the same # of iterations?

for all elasticity tests switching from GMRES->CG yields exactly the same # of iterations. am confused.

Implement MueLu preconditioned solves

The interface already exists, but the implementation is empty. Be sure to pass DOF coordinates to the preconditioner (as this is known to greatly improve preconditioner performance). What should be done for hierarchical DOFs that are not explicitly associated with a coordinate?

Pass coordinates to MueLu preconditioner

This might be a good starting point. This will be slightly complicated by the fact that fields can be mixed order.

elasticity: implement momentum residual evaluator

Balance of linear momentum in absence of body forces.

Consider adding a statistics class somewhere

It's a pain to parse standard output for information like

J(uh) = 1.2984

and

J(eh) ~ 0.003434

Especially because quantities like this are output at different times during the simulation. It would be nice to have some statistics class that can gather information and optionally print a nice data file and optionally print a nice std output. Additionally, it would be cool if this class could optionally accept a value for the exact QoI to compute and output additional statistics (like the exact error and the effectivity of the estimated error) in cases where the exact QoI is known.

poisson: error estimate regression test

check that:
\sum_{i=1}^{n_{vtx}} \eta_i = z^T R(u_fine)

`// namespace ...` to `/* namespace ... */`

All namespaces are currently closed with

// namespace

and all other comments are currently c-style /* */ comments.
This is weird.

poisson: mini-application documentation page

Populate the poisson mini-app doxygen page with meaningful information.

Implement various preconditioner+solver interfaces

Perhaps a few functions like:

solve_ilu_cg(RCP<ParameterList> p, RCP<Matrix> A, RCP<Vector> x, RCP<Vector> b);
solve_ilu_gmres(RCP<ParameterList> p, RCP<Matrix> A, RCP<Vector> x, RCP<Vector> b);
solve_multigrid_cg(RCP<ParameterList> p, RCP<Matrix> A, RCP<Vector> x, RCP<Vector> b);
solve_multigrid_gmres(RCP<ParameterList> p, RCP<Matrix> A, RCP<Vector> x, RCP<Vector> b);

poisson: implement dual solve

The ability to solve problems of the form:
(\nabla z, \nabla w) = (w, j)
in an enriched finite element space.

implement scatter qoi evaluator

scatter dJdu from element-level contributions to the global solution derivative vector.

introduce goal assert macros

Four assert macros should be defined:

GOAL_ALWAYS_ASSERT(cond)
GOAL_ALWAYS_ASSERT_VERBOSE(cond, msg)
GOAL_DEBUG_ASSERT(cond)
GOAL_DEBUG_ASSERT_VERBOSE(cond, msg)

where GOAL_DEBUG_ASSERT* are not called when -DNDEBUG is defined. These macros should replace all instances of c assert in the code.

Error estimation APIs

(1) sum all error contributions from physics->get_e()
(2) sum absolute values of error contributions from physics->get_e() ( an approximate upper bound )
(3) create a APF field from all fields contained in physics->get_e() that is
apf::Field* e_tot = \sum_{i=1}^{nf} abs(physics->get_e()[i])

Implement scatter functional (QoI) evaluator

poisson: verify sum of localized error contributions equals total estimated error

I've already done this, which is super great!
However I mostly hacked it into poisson_solver.cpp.
Maybe I could make this a little cleaner and add a
regression test that checks:
\sum_{i=1}^{n_{vtx}} \eta_i = z^T R(u_fine)

This also indicates that it might be useful to have the following library-wide APIs:
(1) sum all error contributions from physics->get_e()
(2) sum absolute values of error contributions from physics->get_e() ( an approximate upper bound )
I should open a new issue for this....