I am trying to solve the nonlinear Poisson equation of the form
$\frac{d}{d x} \left[ \left( \frac{du}{dx} \right)^2 \right] = f$
where $f=2x$ and $u = \frac{1}{4}$ at $x=0$. I believe the exact solution is $u=\frac{1}{2} x^2 + \frac{1}{4}$. I have the following UFL file:
cell = interval
scalar = FiniteElement("Lagrange", cell, 2)
phi = TestFunction(scalar)
u = Coefficient(scalar)
f = Coefficient(scalar)
F = -(u.dx(0)**2.0)*phi.dx(0)*dx - f*phi*dx
du = TrialFunction(scalar)
J = derivative(F, u, du)
and C++ code:
class Source : public Expression {
public:
Source() : Expression() {}
void eval(Array<double>& values, const Array<double>& coord) const {
const double x = coord[0];
values[0] = 2.0*x;
}
};
class LeftBoundary : public SubDomain {
bool inside(const Array<double>& coord, bool on_boundary) const {
return abs(coord[0])<DOLFIN_EPS && on_boundary;
}
};
int main(int argc, char *argv[]) {
// create mesh
auto mesh = make_shared<UnitIntervalMesh>(100);
// create function spaces
auto V = make_shared<NonlinearPoisson::FunctionSpace>(mesh);
// boundary condition
auto left = make_shared<LeftBoundary>();
auto u0 = make_shared<Constant>(0.25);
DirichletBC leftBC(V, u0, left);
// source
auto u = make_shared<Function>(V); // unknown velocity
auto f = make_shared<Source>();
// create residual form --- (nonlinear) stress balance equations
NonlinearPoisson::LinearForm F(V);
F.u = u;
F.f = f;
// create Jacobian form J = F' (for the nonlinear solver)
NonlinearPoisson::JacobianForm J(V, V);
J.u = u;
// solver parameters
Parameters params("nonlinear_variational_solver");
Parameters newtonParams("newton_solver");
newtonParams.add("relative_tolerance", 1e-6);
params.add(newtonParams);
// solve for the initial velocities
solve(F==0, *u, leftBC, J, params);
}
which complies fine. However, I get a runtime error:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)
Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)
Newton iteration 2: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)
Newton iteration 3: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)
Does anyone know why the solver is producing 'nan's?
