weighted residuals evaluation


106
views
0
3 months ago by
ashoka  
To compute a posteriori error estimators it is needed to evaluate weighted face/element residuals. Before I used:
h_faces = FacetArea(mesh)
h_elems = CellVolume(mesh)

But the problem is that dolfin-adjoint can not derive gradients of these functions.
Is there any effective built-in functions instead?
Thanks!

Community: FEniCS Project
You want to derive gradients from the error estimators? I suggest you use pyadjoint and derive your own Block. If you give more details of your code or a MWE, I can help you write it. More information at http://www.dolfin-adjoint.org/en/latest/documentation/custom_functions.html
written 3 months ago by Miguel  

Yes, I want to derive gradients from error estimator. Actually, I want to implement the topology optimization problem with regularization by a posteriori error estimator. The code looks as follows (it is just the tutorial on linear elasticity + few strings dolfin-adjoint code to derive gradients):

from fenics import *
from fenics_adjoint import *
import pyipopt

# Scaled variables
L = 3; W = 1
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
lambda_ = 1.25
g = gamma

N=10
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(W, L, L), 3*N, N, N)
normals = FacetNormal(mesh)
h_faces = FacetArea(mesh)
h_cells = CellVolume(mesh)

V = VectorFunctionSpace(mesh, 'P', 1) # admissible solutions space
P = FunctionSpace(mesh,'DG', 0) # admissible designs space (piece-wise constant)

# initial design (with respect to which gradient should be computed)
x = interpolate(Constant(1), P)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(x*sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Define cost functional: compliance + a posteriori error
J = assemble(dot(u,f)*dx + h_cells**2*(f + div(x*sigma(u)))**2*dx + h_faces('+')*jump(x*sigma(u),normals)**2*dS)

dJdx = compute_gradient(J, Control(x))

It works well for h_cells and it used to work well for h_faces('+'), but now it gives the following error:

Solving linear variational problem.
Traceback (most recent call last):
  File "t.py", line 58, in <module>
    J = assemble(inner(u,f)*dx + h_cells**2*(f + div(x*sigma(u)))**2*dx + h_faces('+')*jump(x*sigma(u),normals)**2*dS)
  File "/anaconda3/envs/dolfin-adjoint/lib/python3.6/site-packages/fenics_adjoint/assembly.py", line 23, in assemble
    block = AssembleBlock(form)
  File "/anaconda3/envs/dolfin-adjoint/lib/python3.6/site-packages/fenics_adjoint/assembly.py", line 65, in __init__
    self.add_dependency(c.block_variable)
AttributeError: 'FacetArea' object has no attribute 'block_variable'
Abort trap: 6

Thank you in advance!

written 3 months ago by ashoka  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »