# Tag Info

10

I don't think you can avoid using a tolerance for floating-point comparisons. Error due to round-off, discretization, and so on using floating-point numbers is unavoidable. What I typically do to test FEM code I write is: test the stiffness & mass matrices on a single element to make sure I get local element assembly right, compare against a known ...

7

Geoff has already given an excellent overview, but I wanted to provide another real world look on it. In the deal.II project (http://www.dealii.org/) we run some 7,000 tests with every change in the code base, on multiple platforms. The tests Geoff describes are mostly "integration tests", i.e., they run through a significant part of the code base. You need ...

7

You will want to read up on operator splitting methods. In essence, in every "macro time step" you would treat fast processes by doing many "micro time steps" in one half of the algorithm, and then do a single macro time step for the slow processes in the other half. For higher order, you will want to use what's known as "Strang splitting".

6

Off the top of my head, there are a number of things that may be going wrong. You might want to verify the following Are the function values or derivatives you are computing numerically stable? Is the rounding error in their evaluation smooth? You can verify this visually by plotting an extremely small interval around a known root such that the function ...

6

Reformulating (1) and (2), the system reads $$\begin{pmatrix} \mathcal{L}&\\ &\mathcal{L} \end{pmatrix} \begin{pmatrix} W\\ h_2 \end{pmatrix} = \kappa \begin{pmatrix} -\mathcal{L} &\mathcal{M}\\ \mathcal{M} & -\mathcal{L} \end{pmatrix} \begin{pmatrix} W\\ h_2 \end{pmatrix}$$ where the calligraphic symbols are composed of some linear ...

6

The term that you want to search for is multiple timestepping (see, for instance, [1-3]). [1] http://www.cs.unc.edu/Research/nbody/pubs/external/Berne/tuckerman-berne-rossi91.pdf [2] http://www3.nd.edu/~izaguirr/papers/newM3paper.pdf [3] http://arxiv.org/abs/1307.1167

6

There are so many Runge Kutta methods, including Dormand-Prince 45 Cash-Karp 54 Fehlberge (sic) 78 Is there any comparison between them? Well, sure. Here are some traits to compare: Is the method implicit or explicit? (All of your examples are explicit RK methods.) What is the order of convergence? Are there any embedded error ...

6

There are many different ways to do this. One of the standard is a work-precision plot where you plot the amount of time or function calls that it takes in order to achieve a certain level of accuracy. You can find tons of examples at DiffEqBenchmarks.jl. Generally you slide a timestep or adaptivity tolerances along a window and plot all of the (time,error) ...

5

Your question really doesn't admit a simple answer- we need to know more specifics about your problem to provide a useful answer. In general, iterative methods can be faster than direct factorization for large sparse systems of equations if the system is reasonably well conditioned or if it is badly conditioned but you have a good preconditioner or if you ...

5

Daniel Shapero's answer is excellent, but I felt I should add the following: Properly preconditioned iterative methods will almost certainly win here, unless your systems have a very, very special structure. There's a considerable recent literature, unfortunately mostly of the "asymptotic a priori bounds on running time" variety, on solvers for graph ...

5

What do people suggest for the linear system of this type? I know about trilinos, petsc, and sundials, but don't know the other alternatives or have exposure to them. Generally speaking PETSc, Trilinos, and KINSOL (from SUNDIALS) are the best-of-breed when it comes to scalability. From an ecosystem standpoint, PETSc seems most flexible, since it does not ...

5

I am a PETSc developer so take my suggestion with a grain of salt, but I would use PETSc because the problem sizes are large enough that execution overhead should be minimal, you can trivially switch between various sparse and dense solvers (1/2 sparse should be treated as dense, but it might pay off to use a sparse solver for 1/9 at your sizes), a suite of ...

5

If your problems are small, you could try a simple solver with dense data structures and use stack allocation (old Fortran codes). IPOPT uses sparse data structures and indeed, it takes a relatively long time to initialize IPOPT if your problems are small. Old Fortran codes with stack allocation are surprisingly fast if you have small problems. If the ...

5

What is typically done for ODEs is to solve at each timestep using two methods of different order. The discrepancies between the results are then used to estimate the error and then decide whether to decrease or increase the timestep length. You can find more information on the subject here. Two examples of these types of schemes are RKF45 and Dormand-...

5

for 1-D transport problem, under implicit method, we use courant number, a dimensionless number to choose appropriate time steps: $\frac{u\Delta t}{\Delta x}\le C_{max}$ $C_{max}$ should be less than 1, if we make it as 1, we get $\Delta t\le \frac{\Delta x}{u}$ The basic ideal is that, at one time step, a particle placed in the problem domain should ...

5

Since your problem is small, you're probably best off trying fsolve or root. Both of these are interfaces to MINPACK and call HYBRD or HYBRJ. Since calculating a Jacobian matrix for your system shouldn't be hard (either do it by hand, or use your favorite computer algebra system, like SymPy, Sage, Maple, or Mathematica), you should supply a Jacobian matrix. ...

5

I highly recommend the following read: J.R. Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" In short, if the matrix is non-positive definite, there is no guarantee that CG will fail. It might be able to solve it (for some RHSs and certain tolerances), but it is just not supposed to and, probably, converges slower ...

4

Using a line search or trust region method along with Newton's method typically does the trick of finding the zero reliably. Check any advanced numerical methods or optimization book.

4

Can you bound your zero? If you can come up with a good heuristic for the approximate location of the zero and then compute bounds for it over an interval that is not too large, the performance of many methods is improved. Here's some Matlab code I use to perform basic root bracketing for monotonic functions. In terms of interval methods, @Pedro's ...

4

So the problem is that FEniCS doesn't recognize that in 1D, u and grad(u) have the same shape. (You could file an issue about this on the Dolfin issue tracker.) If you replace grad by div, your code works. (The general form of the inviscid Burgers equation is $\mathop{\mathrm{div}}(u^2) = f$.)

4

Which one was faster when you tried using Eigen? You may also consider Trilinos if you're using C++. By "generic Laplacian matrix", do you mean a finite difference / finite element discretization of the Laplace operator on some domain, or do you mean the Laplacian of some graph $G$? Whether or not your solver is overwhelmed as you scale up the number of ...

4

It is quite straightforward to demonstrate this for implicit Euler (aka backward Euler) on a scalar example. Consider the initial value problem $\dot{y}(t) = i \alpha y(t), \ y(0) = 1$ with solution $y(t) = \exp(i \alpha t) = \cos(\alpha t) + i \sin(\alpha t)$. The solution is a harmonic oscillation with amplitude 1 and this amplitude does not change ...

4

Many numerical tips and theoretical explanations can be found in this book from Hairer and Wanner: https://www.springer.com/gp/book/9783540566700 In this book, a strategy is described, which uses a time step such that the relative variation of the solution during the first time step is below a certain threshold if you were using explicit Euler (omitting the ...

4

So so many places you have to rewrite. The whole mesh handling (accessing faces and edges from cells, neighbors from cells, ...). Shape functions. Dealing with the question of how the normal vector of a face when seen from one cell matches that on a neighboring cell. You will also likely encounter that 3d problems are always much larger and that solvers, ...

3

A few additional points I would like to add to other answers. Corner test cases should be part of the regression test suite - ill conditioned problems, ill conditioned - bad aspect ratios, orthotropic and un-isotropic material properties, improperly constrained models. Make sure reaction forces match. Sturm checks for eigenvalue problems especially when ...

3

For a sparse parallel solver, it's your own responsibility to provide a matrix vector product and a suitable preconditioner. The data for the vector itself should fit into main memory in any case. If the matrix has at most a fixed small number of non-zero elements (<20) per column or row, then the same is also true for the matrix itself. In this case, an ...

3

To transform your original program into the form you specified, use the following mappings: First, $D = 2S'S$, and $D$ is positive semidefinite. $b = w$, as you pointed out $d = 0$, as Nico notes $A$ takes the following form: \begin{align} A = \left[\begin{array}{c} I \\ e^{T} \\ -e^{T}\end{array}\right] \end{align} $b_0$ takes the following form: \begin{...

3

If your domain is one dimensional you should use: mesh = UnitIntervalMesh(10) For partial derivative in x direction use: u.dx(0)

3

No, ABAQUS is not a general-purpose PDE solver. It solves the PDE for a number of pre-defined applications such as structural analysis or heat transfer in solids. Of course many physical phenomena are governed by the same equation as heat transfer so you can solve those problems by treating them as heat transfer problems in ABAQUS.

3

You can try using SCS, either the direct or indirect solver. SCS uses first order methods, and hence may be able to solve larger problems than second order solvers such as SDPT3, SeDuMi, MOSEK, etc. On the downside, given that it is a first order solver, it can be very slow - but very slow may be better than not at all. Paper: http://web.stanford.edu/~...

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