Auto adapative on mixed element function space

Avishai_Benyamini · October 19, 2020, 7:29am

Hi,

I am trying to use the auto-adaptive solver to locally refine my mesh. I am following this fenics resource,
https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/auto-adaptive-poisson/python/documentation.html. The only difference in my case is the definition of the Goal functional.

In the demo it is defined with respect to a scalar field

V = FunctionSpace(mesh, “Lagrange”, 1)
…

Define function for the solution

u = Function(V)

Define goal functional (quantity of interest)

M = u*dx

In my case V is a mixed space.

P1 = FiniteElement(“Lagrange”, mesh.ufl_cell(), 1)
R = FiniteElement(“Real”, mesh.ufl_cell(), 0)
element = MixedElement([P1, R])
V = FunctionSpace(mesh,element)

So, M=u*dx won’t work.

I tried:

u = Function(u)
tu = TestFunction(u)

Define goal functional

M = inner(u,tu)*dx

But the AdaptiveLinearVariationalSolver crushes in the end. What is the right way to do this?

The warning:
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.

The error:

RuntimeError Traceback (most recent call last)
in
46 solver = AdaptiveLinearVariationalSolver(problem,M)
47 solver.parameters[“error_control”][“dual_variational_solver”][“linear_solver”] = “cg”
—> 48 solver.solve(tol)
49
50 solver.summary()

~/miniconda3/envs/fenicstest/lib/python3.8/site-packages/dolfin/fem/adaptivesolving.py in solve(self, tol)
75
76 # Call cpp.AdaptiveLinearVariationalSolver.solve directly
—> 77 cpp.adaptivity.AdaptiveLinearVariationalSolver.solve(self, tol)
78
79

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to assemble form.
*** Reason: Expecting a scalar form but rank is 1.
*** Where: This error was encountered inside assemble.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

dokken · October 19, 2020, 9:36am

To me it is not clear what you would like to solve here.
As specified in the demo, M has to be a functional (i.e. returning something in \mathbb{R}.
Assembling with a test function will result in something in \mathbb{R}^M, where M is the dimension of the function space.
How about having the functional M = inner(u,u)*dx?

Note that the updated demo (with latest dolfin syntax can be found at: https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/auto-adaptive-poisson/demo_auto-adaptive-poisson.py.rst

Avishai_Benyamini · October 19, 2020, 3:59pm

Thanks, Dokken!

I initially tried using M=u*dx, but it gave an error as u is a mixed element. So, I tried making it a scalar with a test function. It makes much more sense to use inner(u,u)*dx.

I tried doing that, but I still get an error shown below. FENICS manages to solve the pde without the auto adaptive part.

My code:

Compute solution

w = Function(W)

#Define goal functional
M = inner(w,w)*dx

#Define error tolerance
tol = 1.e-5

#Solve eq with respect to U to minimize the error

problem = LinearVariationalProblem(a,Ls,w)
solver = AdaptiveLinearVariationalSolver(problem,M)
solver.parameters[“error_control”][“dual_variational_solver”][“linear_solver”] = “cg”
solver.solve(tol)

solver.summary()

#Plot
plt.figure(figsize=(20,10))
plot(u.root_node(), title=“Solution on initial mesh”)
plt.figure(figsize=(20,10))
plot(u.leaf_node(), title=“Solution on final mesh”)
plt.show()

Error:

RuntimeError Traceback (most recent call last)
in
46 solver.parameters[“error_control”][“dual_variational_solver”][“linear_solver”] = “cg”
47 solver.parameters[“error_control”][“dual_variational_solver”][“symmetric”] = True
—> 48 solver.solve(tol)
49
50 solver.summary()

~/miniconda3/envs/fenicstest/lib/python3.8/site-packages/dolfin/fem/adaptivesolving.py in solve(self, tol)
75
76 # Call cpp.AdaptiveLinearVariationalSolver.solve directly
—> 77 cpp.adaptivity.AdaptiveLinearVariationalSolver.solve(self, tol)
78
79

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 2 iterations (PETSc reason DIVERGED_INDEFINITE_PC, residual norm ||r|| = 5.708100e+01).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

dokken · October 19, 2020, 4:38pm

As I do not know what problem your are actually solving, I cannot be of much more help. Have you tried to solve the problem without the adaptive meshes?
Does it actually make sense to do the inner product of the full mixed space, or should it be of only the part that is in the vector Lagrange space (not the lagrange multiplier).
Please provide?: A minimal working code example encapsulated with ``` that reproduces the error.

dokken · October 19, 2020, 7:57pm

So your code produces this for me (when changing the tolerance to 1e-6).

Summary of adaptive data:

  Level  |  functional_value  error_estimate  tolerance  num_cells  num_dofs
  --------------------------------------------------------------------------
  0      |          0.104917        0.000328   0.000001        128        82
  1      |          0.093985        0.000081   0.000001        200       122
  2      |          0.088612        0.000014   0.000001        411       236
  3      |          0.085956        0.000004   0.000001        838       462
  4      |          0.084640        0.000001   0.000001       1606       867

Time spent for adaptive solve (in seconds):

  Level  |  solve_primal  estimate_error  compute_indicators  mark_mesh  adapt_mesh    update
  -------------------------------------------------------------------------------------------
  0      |      0.001459        0.003330            0.002701   0.000016    0.000226  0.001664
  1      |      0.000990        0.005357            0.003714   0.000023    0.000439  0.002933
  2      |      0.001508        0.011462            0.006779   0.000048    0.000872  0.005087
  3      |      0.002606        0.022809            0.012928   0.000096    0.001566  0.009134
  4      |      0.004848        0.042813                   0          0           0         0

Isnt this what you expect?

Avishai_Benyamini · October 19, 2020, 8:32pm

That’s what I would expect

For some reason on my system, I get an error for the same code. Is there a reason why this would be system dependent? I am using a mac and running on JupyterLab.

Any thoughts?

---------------------------------------------------------------------------

RuntimeError Traceback (most recent call last)
in
44 solver = AdaptiveLinearVariationalSolver(problem,M)
45 solver.parameters[“error_control”][“dual_variational_solver”][“linear_solver”] = “cg”
—> 46 solver.solve(tol)
47
48 solver.summary()

~/miniconda3/envs/fenicstest/lib/python3.8/site-packages/dolfin/fem/adaptivesolving.py in solve(self, tol)
75
76 # Call cpp.AdaptiveLinearVariationalSolver.solve directly
—> 77 cpp.adaptivity.AdaptiveLinearVariationalSolver.solve(self, tol)
78
79

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 2 iterations (PETSc reason DIVERGED_INDEFINITE_PC, residual norm ||r|| = 9.371859e+01).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

dokken · October 19, 2020, 8:35pm

It shouldn’t be dependent, except if you are using the latest version of PETSc (3.14), which contains quite a few bugs.

Have you tried using docker?

docker run -it -v $(pwd):/home/fenics/shared -w /home/fenics/shared --rm  quay.io/fenicsproject/stable:latest

Avishai_Benyamini · October 19, 2020, 8:39pm

Thanks a lot for helping me debug this.

Actually, now it runs for me with a unit mesh which I probably tried wrongly before.

I am still not able to run this on an arbitrary mesh I am using which I am adding below.

# Create mesh
N=128

##### New Mesh #####
#dimensions
L=4.0
Ws=1.0
a0=4
Wch=0.1
xA=L/2-np.sqrt((Ws-Wch)/(2*a0)).item()
xB=L/2+np.sqrt((Ws-Wch)/(2*a0)).item()
Nch=50 #points on channel
xch=np.linspace(xA,xB,Nch)

#Build domain
#source
domain_vertices=[Point(0,Ws),Point(0,0)]
                                          
# yA=-a*(xA-L/2)**2+(Ws-Wch)/2+0.5

# domain_vertices.append(Point(xA,yA))

#channel down
for x in xch:
y=-a0*(x-L/2)**2+(Ws-Wch)/2
domain_vertices.append(Point(x,y))

#drain
domain_vertices.append(Point(L,0))
domain_vertices.append(Point(L,Ws))

#channel down
for x in np.flip(xch):
y=a0*(x-L/2)**2+(Ws+Wch)/2
domain_vertices.append(Point(x,y))

# domain = Polygon([Point(0,0),Point(L/2-xb,0),Point(L/2,Wsd/2-Wch/2),Point(L/2+xb,0),Point(L,0),Point(L,Wsd),Point(L/2+xb,Wsd),Point(L/2,Wch/2+Wsd/2),Point(L/2-xb,Wsd),Point(0,Wsd)])
domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,N)

plt.figure(figsize=(20,10))
plot(mesh)
plt.show()

dokken · October 20, 2020, 8:26am

I would suggest changing your functional to:

# Define goal functional
M = inner(w[0],w[0])*dx

as you really do not want to include the Lagrange multiplier in the functional.
Then your solver will not diverge

kamilova-maths · March 10, 2021, 3:50pm

Hello,
I am encountering the same exact problem but I have not been able to fix it with the discussion here, even when changing my mesh for a unit mesh and running it using docker.
Did you manage to implement it on an arbitrary mesh? Could I see the working code?
Thanks so much!

Topic		Replies	Views
Mesh adaptivity in fenics	2	1649	August 29, 2023
Problem with Adaptive LinearVariationalSolver Errors	1	625	May 19, 2021
Error on the goal function(?) in the Adaptive linear solver	1	416	March 11, 2020
Solving adaptive nonlinear equation in Vector Space variational formulation	11	1829	February 20, 2021
How to return the data of the adaptive solver? I/O	0	390	September 25, 2019

Auto adapative on mixed element function space

Define function for the solution

Define goal functional (quantity of interest)

Define goal functional

The error:

Compute solution

Error:

Related topics