Hot answers tagged

11

Julia is built in such a way that you will never see a full PETSc-like library, and that's on purpose. PETSc is not a single thing: it is an HPC library with some utility functions, linear solvers, nonlinear solvers, ... a whole can of soup that works within its own world but not outside of it. Julia's package ecosystem is build on generics and composability....


10

I've been developing a lightweight finite element library in Python 2.7 harnessing the power of NumPy arrays and SciPy sparse matrices. The general idea is that given a mesh and a finite element, you have more-or-less one-to-one correspondence between the bilinear form and a (sparse) matrix. The user can then use the resulting matrix as he or she sees fit. ...


8

As one of the library's authors, I would of course love for deal.II to come out on top with this comparison. But I suspect it may not, and the answer lies in a factor you omit from your comparison: how long it actually takes to implement your code. Few people in academia with the skills to implement a FEM code from scratch spend more time solving PDEs than ...


6

I would recommend looking at a time-dependent example because PETSc can provide a lot more diagnostics if you formulate at that level. For example, you could use a Rosenbrock method, an additive Runge-Kutta IMEX method, or others. Some involve Newton iteration, but that is not required and is not always the best approach. To use those methods, you can ...


5

Some thoughts from someone who has worked a fair amount in compiled languages, and has done a tiny bit of FVM: Typically, if you have experience programming in C, you sketch out a high-level description (pseudocode) of what you would like to do. Then you look for libraries that might implement the data structures and capabilities you need for your high-...


5

You say that you want an MPI version. Then you need to study the literature, as the distributed memory variant of matrix-matrix product are not a simple parallellization of the sequential version. The Cannon algorithm is pretty cute if you're on a square processor grid. In each step you rotate the input matrix rows and columns, so that in the end each ...


5

Your problem is too small. You have to consider that to get good efficiency, each processor has to have enough work to offset the cost of communication. In other words, there is a threshold how many degrees of freedom there have to be below which efficiency deteriorates. I don't know where that threshold is for MUMPS. But, to give just one example, for ...


4

Unfortunately there's no tool for this. You can run each on a variety of input sizes to establish the computational complexity they appear to have, i.e. the $f$ in the $O(f(n))$ that characterizes each code. This can point you towards what underlying algorithms each is using and verify or not that something that should be $O(n)$ is actually implemented to ...


4

If you're using iterative methods with assembled matrices, just buy DDR channels. Don't pay attention to number of cores when loooking at the spec sheet. Within the same class of processors (e.g., a recent generation of Xeon), the achievable memory bandwidth will be proportional to the quoted peak bandwidth. Note that for very small problem sizes, you might ...


4

I like this formulation for implementing Dirichlet boundary conditions in cases where elimination of boundary dofs is not convenient. If you apply effectively the same procedure to $M$, the rows and columns corresponding to Dirichlet boundary dofs will be rows and columns of the identity. This will give you new eigenvalues equal to 1 (or $\alpha$ in the ...


4

Continuous weak form Though I think the weak form is more fundamental, suppose we start with the strong divergence-form representation for a (first- or) second-order quasilinear PDE: find the $m$-component solution $u \in R^m(\Omega)$ $$ -\nabla \cdot f_1(u,\nabla u) + f_0(u,\nabla u) = 0 $$ on the domain $\Omega \subset R^d$ with (for simplicity) ...


4

It is indeed possible to write an obstacle problem for an advection-diffusion equation as a variational inequality: If $a(u,v)$ is the bilinear form corresponding to your advection-diffusion equation, the corresponding obstacle problem (in your specific case) is finding $u\geq 0$ such that $$ a(u,v-u) \geq (f,v-u) \qquad \text{for all } v\geq 0.$$ In ...


4

I think you have some confusion. PETSc is not in the same league as Fenics, Libmesh, Moose etc. In fact, all of these (heavyweight) packages use PETSc for linear algebra. IMHO PETSc is as lightweight as you can get. It just requires C/Fortran compilers and Python (used only for configuration) and you can build the library in under 5 minutes on your laptop. ...


4

I can recommend nutils. nutils meets at least a few your "light-weight" requirements. it is pure python and easy to install since it only depends on standard Python libraries numpy, scipy, and matplotlib and, thus, it is well suited for interoperations. At least the developers claim that "Exposed objects are of native python type or allow for easy ...


4

Wow! A question that I can answer! I have been using PETSc for the past year and a half to solve the Navier-Stokes equations (with some hard-coded MPI). The best way to learn PETSc is to (1) read the manual so you generally know what's in PETSc, (2) decide roughly which part of PETSc you need, and then (3) begin experimenting with simple examples to ...


4

Such problems can be solved using an $LDL^T$ factorization (similar in memory and time cost to Cholesky). Your matrix is not very sparse so treating it as such may have limited benefit. I would recommend comparing sparse direct solvers (such as MUMPS and Umfpack) to the dense factorizations you've used. Complex-symmetric matrices are often difficult to ...


4

Use of Lagrange multipliers produces a saddle-point problem, $$ \begin{pmatrix} A & B^T \\ B & 0 \end{pmatrix} \begin{pmatrix} u \\ \lambda \end{pmatrix} = \begin{pmatrix} b \\ 0 \end{pmatrix} $$ As you've noticed, many preconditioners break down for this sort of system. One can use direct solvers that support pivoting, but if you want iterative ...


3

As mentioned in Geoff Oxberry's more complete answer, it should be noted that PETSc includes TSSUNDIALS, an interface to SUNDIALS. If you configure PETSc with the --download-sundials option (see python2 ./configure --help | grep -A 2 sundials for other related options, such as using an existing SUNDIALS library), then you can use at least some of the ...


3

In general, you can do more with PETSc. SUNDIALS is a collection of ODE solvers (in CVODE, Adams-Bashforth and BDF methods; in ARKODE, ARKIMEX methods) and DAE solvers (IDAS implements a BDF method) with sensitivity capabilities (the CVODES and IDAS variants), and a nonlinear solver (KINSOL). There are a few Krylov solvers in there (at least GMRES and CG) ...


3

The Teuchos Unit test harness in Trilinos natively supports unit tests that use MPI. Things like controlling output from multiple processes and aggregating pass/fail over all processes is automatic. Take a look: http://trilinos.org/docs/dev/packages/teuchos/doc/html/group__Teuchos__UnitTest__grp.html


3

ADI is not a very good parallel algorithm. You should seriously consider formulating the problem in 3D and solving with multigrid. You could get a start with src/ksp/ksp/examples/tutorials/ex45.c. If you insist on using ADI, you should seriously consider allocating the three matrices separately. It's more memory, but then you won't have to reassemble on ...


3

As of writing this answer (October 2014), Trilinos does not have this capability throughout its code base. There are two packages (Kokkos and Tpetra) that will provide this functionality to Trilinos, but neither are used widely throughout the remainder of the library. The Trilinos plan (as far as I understand it) foresees new packages (and some revamped old ...


3

Since nobody is stepping up to answer this question, let me point to a paper that contains an example of why one would want to do this: http://www.math.tamu.edu/~bangerth/publications.html#publications/reviewed-34 (publication 34, Kronbichler et al.). In fact, there are multiple examples in this paper: We precondition $\begin{pmatrix}A &B \\ B^T & ...


3

Krylov solvers have no problem converging if the system is consistent, i.e., if the right-hand side is in the image of the system matrix $A$; see, e.g., Iterative Krylov Methods for Large Linear Systems by van der Vorst. PETSc's CG solver should work. You can even use a multigrid preconditioner if you make sure that the coarse solver preserves the null ...


3

If you're looking for a book and are happy with finite element discretizations, you can check out Elman, Silvester, and Wather, "Finite Elements and Fast Iterative Solvers". If you prefer finite volume or finite difference methods, Wesseling's CFD book has a section on multigrid that will at least give you the context for multigrid-specific material. If ...


3

There are many methods for solving linear systems in parallel, but fundamentally what works and what doesn't depends crucially on the properties of the linear system you are trying to solve. Among those, the Jacobi preconditioner is about the worst you can choose. I've summarized what I know about this topic in lectures 34-38 and 41.75 here: http://www.math....


3

I don't think PETSc supports it. PETSc really thinks in parallel, so it converts at most between distributed matrices and sequential matrices through submatrix taking operations. I would MatGetArray and send that with an MPI call.


2

If you have an idea about how to statically partition your problem into subgroups, you could try partitioning the MPI processes into these subgroups, and then creating an MPI communicator for each subgroup. Most PETSc & SLEPc object creation routines take a communicator as an argument, so instead of using MPI_COMM_WORLD, you could use the ...


2

What if one Vec was created by VecCreateSeqWithArray and the other not? How would the user know which Vec holds the memory they allocated? There is a dirty backdoor in VecPlaceArray, but using it is typically bad design.


2

It depends on the object representation of your 3D array. If it's a PETSc Mat, you want to look at MatView; if it's a Vec, you want to look at VecView. The links are to PETSc documentation, and each of the links has a list of relevant examples. They may not be exactly what you're looking for, but they should be close; in general, you're looking for something ...


Only top voted, non community-wiki answers of a minimum length are eligible