[KTG]

Hello,


I am trying to solve the emden-fowler equation: T*laplace(T)=x by writing an openfoam solver (I only care about the 1D problem, so I am taking T as a scalar). I have some code that can solve div(T)=x using fvOptions - that part is not hard. However, when I add in multiplication by T, it complies, but the solver crashes. I am not sure what is wrong:

                fvScalarMatrix TEqn
                (

                       T*fvm::laplacian(DT, T) == correctDim*fvOptions(T)
                );

                fvOptions.constrain(TEqn);
                TEqn.solve();
                fvOptions.correct(T);

or if the solver is actually the problem. When I remove the fvOptions forcing term, the solver will run, but produce incorrect results. For example, if I specify dirichlet 1 boundaries on either side, the solution is f(x) = 0 with kinked boundaries, but should be f(x)=1.


Any ideas? I was thinking maybe I am not treating the T*laplace(T) term correctly?


Thanks

[KTG]

So here is a little update: It looks like I was actually crashing the matrix solver. So I think I am solving the correct equation, but now I am getting a singular matrix. I am solving:

                fvScalarMatrix TEqn
                (

                       T* (fvm::laplacian(DT, T)) == -2*correctDim*fvOptions(T)
                );

                fvOptions.constrain(TEqn);
                TEqn.solve();
                fvOptions.correct(T);

with x in [.2 1], fvOptions is coded to be just the spatial values of x, and the left boundary is zerogradient while the right is fixedvalue 0. Its tempting to say that its probably a boundary condition issue, but I know this specific problem to have a published power series solution... Anyone have any ideas? ANY ideas... Like is there a way to print the matrix problem from the solver so I can manually check it for problems? Given how simple the problem is I probably could have written in house code for this by now (which is not to say that would work)...



Thanks