[benk]

Hi, I've been trying to figure this out for some time now. I have experience using Comsol and also OpenFOAM.

Comsol seems to be able to solve the following type of 1D coupled problem on a domain x=[0:1] quite easily:

PDE #1: div( D1 * grad(c1) ) = -(c2-c1)
Boundary conditions:
at x=0 D1*d(c1)/dx = -1 (flux boundary condition)
at x=1 d(c1)/dx = 0 (zero gradient)

PDE #2: div( D2 * grad(c2) ) = (c2-c1)
Boundary conditions:
at x=0 d(c2)/dx = 0 (zero gradient)
at x=1 c2 = 0 (fixed value)

(these are the parameter values that I used: D1=1000, D2=1e-2)


At first glance this might seem simple enough (perhaps it is), but PDE #1 has two flux boundary conditions and so its solution could basically lie anywhere as long as the boundary conditions are met. Yet, when solving this in comsol, I get a solution that seems to make sense. Note that when solving PDE #1 without any coupling (ie. setting the source term to a constant value), the solution is very different.

So my question is:
How does Comsol (or any other pde solver) figure out what a proper reference value for c1 should be in this case? I'm not providing any reference value yet there is a solution to the above hypothetical model.

Are there any good reference material that I can read concerning the solution of coupled pde's?

[Rami]

Hi Benk,

If your PDEs are decoupled (e.g., use RHS=0 for the first, as you suggested) then the solution of PDE1 is indeed "floating" (i.e., determined only up to a constant) because of the two derivative BCs. However, since the PDEs are coupled, and since PDE2 has a Dirichlet BC, it fixes the solution for both.

I am not familiar with COMSOL, so I can only guess that it gives you just one of the infinite soultions for the uncoupled PDE1, which probably is correct. If there is a way to fix one point within the domain in this case, the solution should be unique.

Merry Chrismas and a Happy New Year,
Rami

[benk]

Thanks for your reply.

I'm wondering how this is actually implement in other codes. For example, in OpenFOAM, we can set a reference point for a floating solution. So I suppose a first estimate of this reference point for PDE 1 would be 0 (the fixed value of pde 2). Does this make sense?

But then, setting the reference point to 0 won't give a converged solution, you'll have to essentially search for a better reference point by moving the reference point up or down until you find a converged solution. Are there any established "searching algorithms" for this sort of thing?

[Rami]

Hello Ben,

I can't tell you how other codes are solving, so I'll stay on the principle level.

Assuming your problem is well-defined (see my former reply), it can be solved directly (if linear) or iteratively. For the direct solution, there is no need for an initial guess (or reference value). For the iterative solution, you usually prescribe an initial guess (so the algorithm does not need to supply one for you) and then the iteration is taking care of updating the solution until convergence. If the algorithm is robust and the problem is well-behaved, it will not be much affected by your initial guess, except for doing more iterations for a poor guess. Otherwise, it is your responsibility to try to find a good initial fields. To this end, you may linearize the equations, use degenerate forms or use your gained experience.

I hope I answered your question,
Rami