19
There are two issues that you are likely to be encountering.
Ill-conditioning
First, the problem is ill-conditioned, but if you only provide a residual, Newton-Krylov is throwing away half your significant digits by finite differencing the residual to get the action of the Jacobian:
$$ J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$$
If you ...
10
Let me first answer all the questions:
What is the theoretical convergence rate for an FFT Poison solver?
The theoretical convergence is exponential as long as the solution is sufficiently smooth.
How fast should this energy converge?
The Hartree energy $E_H$ should converge exponentially for a sufficiently smooth solution. If the solution is less ...
9
Just as an aside, your github documentation is fantastic.
This is just a guess from DG methods, which can have similar issues if numerical fluxes aren't chosen carefully (I figure FV methods are a subset of DG methods). If you're using interpolation from cell centers to define your fluxes, then this should be equivalent to using the average as a numerical ...
8
Let's examine the one-dimensional three-point stencil case in detail,
because I think it's important to be clear just how this behaviour
arises, and what it means to set a point to a certain value in a
finite-difference grid when the underlying function is
discontinuous. The equation will be
$$ u''(x) = \rho(x). $$
Instead of using the interval $[-1,1]$ with ...
8
When looking at the solution of your system, you will find that almost all entries of $x$ are nonzero although the right-hand side is "sparse". Hence, whatever algorithm you use, it'll have to visit each and every entry at least once, so one wouldn't expect that you can save a lot of time using the sparsity of $b$.
Right-hand sides where you can save a lot ...
8
It is standard procedure for general case like that. You can enforce condition for electrode by Lagrange multipliers,
$$
L(u,\lambda) = \int_\Omega \sigma u_{,i} u_{,i} \, \textrm{d}V - \int_{\partial\Omega_\sigma} u n_i j_i \textrm{d}S + \lambda\left(\int_E \sigma u_{,n}\, \textrm{d}S - I \right)
$$
that for discretised problem is
$$
\left[
\begin{array}{...
7
I think you are on the right way. If you correct your errors, it will look very similar to http://www.math.toronto.edu/mpugh/Teaching/Mat1062/notes2.pdf.
7
For your example equation, taking the average approach, the local consistency error
$$
\frac{1}{h^2}[u(x-h) - 2 u(x) + u(x+h)]-f(\frac{1}{2}[u(x-h)+u(x+h)]) = \frac{1}{2}f_uu_{xx}h^2 + hot.
$$
will be of order $2$ (instead of order $3$). ($hot.$ means higher order terms)
Therefore, if your overall approximation is of order $1$, e.g. if you use upwind ...
6
Of course you can discretize the two equations of your system with two different methods. The challenge will simply be when the solution of one equation enters that of another. At that point, you will have to decide what the finite volume solution should be at a finite difference point should be, or how to integrate the finite difference solution (which is ...
6
For a 6x6 grid, those are about the error differences I would expect from two different methods. You have to realize that a 6x6 grid is a very coarse grid, even for a simple problem like yours. As long as you see the two solutions converge towards each other as you refine your grid, there is likely no implementation error. Finite-difference has no general ...
6
Preliminary remarks:
The natural jump condition is on the normal derivative, not on the full gradient.
The fact that the solution is continuous across the interface is part of the formulation of the problem. This is not only due to the finite element formulation.
Formulation
The standard formulation simply enforces it with... no extra term.
Start by ...
6
You are solving Poisson's equation:
$$\nabla^2 \Phi = 4\pi G \rho.$$
Notice that if $\Phi$ is a solution, then so is $\Phi+C$ for any constant $C$. Furthermore, $C$ will have no effect on your simulation since the forces depend only on derivatives of $\Phi$.
The divide-by-zero issue is just a manifestation of this degeneracy. In the Fourier-...
5
Your formulation of $A$ assumes that $u_0$ and $u_{n+1}$ are zero, which is correct. However, you are then including your boundaries at $u_1$ and $u_n$. Your exact answer satisfies these later BC's but not the imposed BC's in $A$. These discontinuities floating around cause you to lose your higher order accuracy. To get second order convergence you only ...
5
I think this happened because $L$ is a finite-difference operator that approximates your equation with some error (it is $\mathcal O(dx^2)$). To get small error your should increase nx and ny. If you set in your example nx = ny = 101 ($dx = 0.01$) you will get error about 1e-4.
5
The answers so far seem to me to be completely off target. The boundary conditions you have are consistent with the r.h.s. of the equation and the matrix $L$, so they are not the problem.
Instead, there are two things wrong here. First, the usual finite second-order difference approximation to $u''$ is
$$ u''(x_0) = h^{-2}\big(u(x_{-1}) - 2u(x_0) + u(x_1)\...
5
If you have simple, grid-aligned interior objects or accuracy is not crucial then stick with the method you have described. If you need to accurately represent arbitrary boundary shapes then you're probably best off moving to a (more complicated) finite element approach on an unstructured grid.
You have described the simplest approach to solving this ...
4
Your "straight line" approach would mean that $\phi_1-\frac{\phi_2-\phi_1}{x_2-x_1}(x_1-x_0)=0$, where $x_0$ is the location of the boundary and $x_i$ are the locations of the places where you define the $\phi_i$. This means that you can eliminate $\phi_1$ in terms of $\phi_2$, in the same way as in your first approach you eliminated $\phi_1$ right away by ...
4
Is the sign accurate? If so, you may have an issue since your var form should be
$(\nabla u, \nabla v) - \langle n\cdot \nabla u, v\rangle_{\Gamma_{\rm rest}} = 0$
substituting the Robin condition in gives
$(\nabla u, \nabla v) - \alpha\langle u, v\rangle_{\Gamma_{\rm rest}} = -\langle 1, v\rangle_{\Gamma_{\rm rest}}$
which can mess with your ...
4
It's not a question whether the condition holds or not. You choose $R$ so that it's a multiple of $I^T$ because you want the condition to hold. You want the condition to hold because that's the only way you can ensure that the multigrid operator is symmetric and amenable to, say, being a preconditioner for CG.
4
So Is it neccesary to use Poisson - Boltzmann equation if I only need to build electrostatic potential from a PQR file?
No. You can use Poisson.
Since you know the positions of each point charge, you know the charge distribution $\rho$, which is a sum of delta functions.
You can thus solve numerically the Poisson's equation that links the charge ...
4
First, some notation for clarity. We're talking about solving $u_{xx}=f(x)$. So $u$ is the solution and $f$ is the right hand side.
What makes you think either of your solutions is or is not correct?
To check your results, just compute two derivatives of your solution (using finite differences or FFT, either one) and compare with $f$. In your example, $f$ ...
4
As noted by others, this PDE is not well-posed in the Strong Form. However, you can take some function which asymtopically approximates the Dirac Delta and use that (where the error in the approximation is below truncation error). For example, you can use
$$ \eta_\epsilon = \epsilon^{-n} \eta(x/\epsilon) $$
where
$$ \eta(x) = \exp(-1/(1-|x|^2)) \chi_{|x|&...
4
For solving systems you shouldn't be comparing the results between each iteration but rather computing the residual.
If you consider the matrix representation to be in the standard form:
$ A\cdot y = b $
Then you can define the residual at some iteration ($y_{i}$)
$ r =b - A\cdot y_i $
Then you 'just' need to determine a stopping method based on the ...
4
The finite difference matrix for the Poisson equation is symmetric and positive definite.
So the preconditioned conjugate gradient algorithm is the iterative solver of choice for
this problem.
The choice of preconditioner has a big effect on the convergence of the method. Incomplete
Cholesky factorization is known to work well for this problem.
The ...
4
Like @likask points out in the other answer, using Lagrange multipliers is one way to do this. The other is to recognize that if you have constraints of the form $\sum_i a_i U_i=0$, then you can select one degree of freedom (typically the one with the largest $a_i$, to maintain numerical stability and rewrite the constraint as
$$
U_\ast = -\sum_{i, i\neq \...
4
$h$ is a measure of the mesh size. In the example, they are using rectangular elements. For which a commonly used measure of the mesh size is the length of the largest diagonal.
Looking at the table above the sentence you quoted:
$(M_5-M_6)/(M_4-M_5)=(0.14052586 −0.14056422)/(0.14037251−0.14052586)\approx 1/4 = (h_6/h_5)^2$
where $M_i$ is the mean at $i$-th ...
3
My suggestion is to just write your own code using PETSc instead of using an existing FE library. Parallel assembly/solve is the most complicated part of an FE code and PETSc takes care of it. The rest of it is simple anyway, maybe a few hundred lines of straightforward C/Fortran. Besides PETSc now has DMPLEX for managing meshes (though it is still in ...
3
The basic equation of electrostatic is
$$
\mathbf{E}=-\nabla\varphi
$$
not
$$
\frac{\partial\varphi}{\partial x}=-\mathbf{E}+C
$$
It follows from a simple principle that the work electric field does when moving a charge from point A to point B does not depend on the path between A and B, otherwise, we would have a perpetuum mobile. Potential is by definition ...
3
Your matrix is correct, but your boundary conditions are not being enforced on your RHS vector.
For example, consider the row corresponding to an interior corner point that will use the two distinct boundary values. The RHS vector of this row should be
$$-1.0 - 1.5$$
because the corners require two values to be subtracted to the RHS. The same is ...
3
In the finite element method, you need to compute matrix entries (here in the case of the Laplace equation)
$$
A_{ij} = \sum_K \int_K \nabla\varphi_i \cdot \nabla \varphi_j
$$
and right hand side entries
$$
F_{i} = \sum_K \int_K \varphi_i(x) f(x)
$$
where $K$ are the cells of the mesh. The former can sometimes be computed analytically, but depending on ...
Only top voted, non community-wiki answers of a minimum length are eligible
