8

Continuous finite elements Typically, if $A$ is your finite element discretization on the finest mesh, $A_i = R_i * A * R_i^T$. So, for $i=0$, $A_0$ corresponds to the finite element discretization of your problem on the coarsest mesh. If $R_0$ and $R_0^T$ are algebraically defined, you can use the "Galerkin" triple product to form $A_0=R_0 * A * R_0^T$. ...


7

If you are not using AMR and do not want to scale beyond 1K-4K cores then simply do this. Rank 0 reads the entire mesh and partitions it using METIS/Scotch etc. (Note: This is a serial operation). Rank 0 broadcasts the element/node partitioning info to all other ranks and frees the memory (used to store the mesh) All ranks read the nodes/elements they own (...


6

Short answer: yes, you have to use something different, e.g. a Neumann-Neumann method. A good reference is Widlund's book. Non-overlapping methods are based on the principle that, if $u$ solves the Poisson equation in a domain $\Omega$ which has been partitioned into two domains $\Omega_1$, $\Omega_2$, at the shared boundary $\Gamma$ of both domains, $u$ and ...


4

The classification of domain decomposition methods into non-overlapping and overlapping methods can be refined. The non-overlapping methods have two groups: single trace and multiple trace methods, where `trace' is sometimes the primal variables (aka Dirichlet traces) and sometimes the dual variables (aka Lagrange multipliers) on the interfaces. In the ...


3

For your first question, constructing the adjacency graph of the "partitions" (what you call "cell groups"): Let's say you have an array $p_K$ in which you store for each cell $K$ which partition $p$ it belongs to. Also assume that you have a (sparse) array $a_{KL}$ whose entries are true if cells $K$ and $L$ are neighbors ("adjacent&...


3

If different nodes have different costs, for example because different rows of your matrix have different numbers of nonzero entries, then you need to attach weights to each node of your graph. Graph partitioning algorithms such as METIS allow you to do this, creating partitions where it is not the number of nodes that are about equal between partitions, but ...


3

$c$ depending on time is not the issue. You will use an RK scheme which takes care of this. The issue is $c$ is discontinuous in $x$. I recommend SBP-SAT schemes for this. (1) Derive an energy equation at PDE level. (2) Search literature for SBP-SAT schemes which enforce interface conditions via SAT penalty terms, which are designed to mimic the energy ...


3

Don't use domain decomposition methods. They're from the 1990s, but we have much better ways of preconditioning problems today. All of them work on the global problem, rather than ones on subdomains. One example are algebraic multigrid methods.


3

Partitioning by elements leads to minor simplifications in assembly at the expense of poor load balancing in the solver, inability to use many popular solution methods, and more communication to apply the operator. The most common reason to use an element partition is simply because the code to handle overlap or to communicate contributions into off-process ...


2

This may not come as a surprise to you given that I develop deal.II, but here's my perspective: When I talk to students, I typically tell them to develop their own prototype in the beginning so they can see how it's done. But then, once they've got something small running, I make them use a library that allows them to go so much further because they don't ...


2

Both are valid techniques with their own tradeoffs. PETSc, in particular, has a history as a linear algebra library, and given that, it makes sense for it to assign rows distinctly to processors. For finite element methods, which have element matrix formation, and often some sort of matrix assembly, it often makes sense to assign elements to uniquely to ...


2

Let us consider the abstract linear problem $$ \mathcal{A} x = b \,, $$ where $\mathcal{A}$ is a linear operator and $x$ and $b$ some functions on a certain domain. To answer your question let me first discuss, what a preconditioner is supposed to do. A preconditioner is choosen such that $\mathcal{B} \approx \mathcal{A}^{-1}$, implying that $\mathcal{B} \...


2

Here is a brute force solution that would work no matter what is the discontinuity and nonlinearity in $c(x,t)$. Write your PDE as a system of two: $ \dot{y}=z\\ \dot{z}=c^2(x,t) y_{xx} $ Now, discretize it on a uniform spatial grid in x: $ \vec{x}= [x_0, x_1,..., x_{n-1}] \\ \vec{y}= [y_0, y_1,..., y_{n-1}] \\ \vec{z}= [z_0, z_1,..., z_{n-1}] \\ $ Now the ...


2

This is always going to be a somewhat costly calculation because the inverse of any "interesting" sparse matrix is generally dense and therefore so is $\mathbf M$. That said, there are smarter techniques than the obvious one you mention, that are asymptotically no slower than the factorization of $\mathbf Q$ itself, as long as $\mathbf B$ has no more rows ...


2

In general, the appropriate preconditioners for elliptic problems such as yours are multigrid methods. In this 1d case, however, the simplest discretizations lead to tri-diagonal matrices and in that case the Thomas algorithm can be used to solve the problem directly without too much trouble. So you don't even need a preconditioner if you use a three-point ...


1

I found the ParMetis have what I want and easy to use.


1

It's a misunderstanding that you need two different meshes: The proper way to see things is that you are using the same mesh, but different polynomial spaces for the two variables. For example, for the Stokes equation, you'd have quadratic polynomials for the velocity $\mathbf u$ and linear polynomials for the pressure $p$. Appropriate parallelization ...


1

Here's the flow of logic, hopefully in a slightly more readable form: We write the transmission conditions with as-yet-unknown linear operators $\mathcal{S}_1$ and $\mathcal{S}_2$ which operate on the fields $u_1$ and $u_2$ at the interface. Perform a Fourier transform in the direction $y$ which is tangential to the interface. $k$ is the spatial frequency ...


1

The block SPD-ness is preserved because any partitioning does not change the SPD property of the matrix as it would not alter the eigenvalues as partitioning would involve row and column swaps which would only re-arrange the matrix but not change the fact that the eigenvalues of its block be greater than 0.


1

As I understand, this is a fairly popular approach. Direct solvers are usually more efficient than iterative solvers for < 100,000 unknowns, so you partition the problem into subproblems of roughly this size, use direct solvers on each subproblem, and combine them globally with some method like alternating Schwarz. It's also common use the domain ...


1

For your example, it is the ratio of the overlap size and the subdomain size that matters. With coarse grid, the condition number scales like H/d with H the subdomain size and d the overlap size. It does depend on the geometry of the subdomains but it seems unaffected from high-order finite elements because this is a property holding already at the ...


1

I'm not sure you need a Steklov-Poincare operator for overlapping domain decomposition methods. The Dirichlet-to-Neumann map (i.e. DtN or Steklov-Poincare operator) is useful in non-overlapping domain-decomposition because Neumann (i.e. flux) data is used in the transmission conditions between subdomains. A discrete DtN map just does this by solving the ...


1

The problem with elliptic operators is that the solution at any specific point depends on all the domain. The overlap allows to 'mix' the approximate solution between subdomains and therefore accelerates the solution of the problem. http://en.wikipedia.org/wiki/Domain_decomposition_methods


1

It looks like the answer to the question is no; however, Japhet, Maday and Nataf have come up with a way to formulate Robin transmission conditions in tandem with a Neumann-Neumann Lagrange multiplier approach. The approach appears to be: initialize $u_0, \lambda_0$ by solving the standard Neumann-Neumann problem with Lagrange multipliers. $\lambda$ ...


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