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 ...
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
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
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
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
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
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
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
One of the standard ways to handle singular problems like Poisson with Neumann boundary conditions which lead to a sinuglar problem $Ax=b$ is to make the obtained solution orthogonal to the kernel (which is for simple domains only the constant vector) if you use something like a Krylov subspace method.
For instance, you can solve your system with conjugated ...
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
Previous comments gave you good suggestions, I try to add some more.
Firstly, your example 2x2 really does not correspond to zero Neumann boundary conditions. In fact, one can show that your choice leads to zero Dirichlet boundary conditions.
If you domain is a rectangle, then implementing the zero Neumann condition can be done easily. You have to include ...
3
Your difference stencil is wrong for nodes at the boundary (which in your 2x2 case happens to be every row). A simple way to see this is to realize that if you applied the Laplace operator to a constant function, the result is zero. Similarly, if you applied the Laplace matrix to a vector where all components have the same value (the discretized version of a ...
3
Assuming you choose stable pairings of elements, the main consideration is which variable you are more interested in. For example, if you choose the $H(div),L^2$ formulation, then you have the choice of piecewise constants for the pressure and a BDM element for the velocity. This yields only first order accuracy for the pressure in the $L^2$ norm -- in other ...
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 stability-analysis of FVM discretizations for elliptic problems with Dirichlet BC, a central assumption is that the inner cells, where you state the PDE, have no intersection with the boundary, i.e.
$$
\overline \Omega_i \cap \Gamma_D = 0 \quad \quad (*)
$$
if seen as a set in $\mathbb R^{n-1}$ if your domain $\Omega \subset \mathbb R^n$, cf., e.g., ...
Only top voted, non community-wiki answers of a minimum length are eligible
