Try to evaluate the sensitivity value but couldn't solve the correct value

Pasuta · April 11, 2024, 4:31pm

I am now solving a topology optimization problem that requires solving for sensitivity by substituting concentration, eps as the design variable and lamda into this equation
.

dJde = (k0*gamma*((1-eps.x.array)**(eta-1))*c.x.array) + (D0*eta*(eps.x.array**(eta-1)))*(grad_c[0]*grad_lamda[0]+grad_c[1]*grad_lamda[1])+(k0*gamma*((1-eps.x.array)**(eta-1))*c.x.array*lamda.x.array)

But the answer results are as below. which seems impossible to calculate

import dolfinx 
import numpy as np
import dolfinx.nls.petsc
import ufl
import dolfinx.fem.petsc
import pyvista
from dolfinx.mesh import CellType, create_rectangle
from mpi4py import MPI
from dolfinx import fem, nls, io, default_scalar_type
from ufl import Measure, inner, grad, div, TestFunction, TrialFunction, adjoint, derivative, SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from petsc4py import PETSc
from dolfinx.plot import vtk_mesh


#Parameter
D0 = 0.1;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;


# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [100, 100], CellType.triangle)
V = fem.FunctionSpace(mesh, ("CG", 1))
c = fem.Function(V)
phi = TestFunction(V)

# The eps is the design variable
eps = fem.Function(V)
eps.x.array[:] = 0.5

#boundary_markers = FacetFunction('size_t',mesh)
dx = Measure("dx") ##integration over different parts of the domain
x = SpatialCoordinate(mesh)
gu = Constant(mesh, default_scalar_type(0)) ##Neauman = 0

# Residual of the variational form of Poisson problem
R = inner(D0*eps**eta*grad(c), grad(phi))*dx - inner(k0*(1-eps)**gamma*c,phi)*dx - inner(gu,phi) *ds     # " - inner(gu,phi) *ds " is boundary neauman 

##Dirichlet 
def on_boundary(x):
    return np.isclose(x[0], 0)
dofs_L = fem.locate_dofs_geometrical(V , on_boundary)
bc_L = fem.dirichletbc(default_scalar_type(1), dofs_L, V) ##Dirichlet = 1
bcs = [bc_L]

#Forming and solving the linear system
problem_u = fem.petsc.NonlinearProblem(R, c, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD,problem_u)
solver.rtol = 1e-16
solver.solve(c)

#how big is the solution
J = -ufl.inner((k0*(1-eps)**gamma),c)*dx ##Objective function
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)

#solve the adjoint vector
lhs = ufl.adjoint(ufl.derivative(R,c))
rhs = -ufl.derivative(J,c) 


def on_boundary_1(x):
    return np.isclose(x[0], 0)
dofs_L_adj = fem.locate_dofs_geometrical(V , on_boundary_1)
bc_L_adj = fem.dirichletbc(default_scalar_type(0), dofs_L_adj, V) ##Dirichlet = 1
bcs_adjoint = [bc_L_adj]

#Adjoint problem
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs_adjoint, petsc_options={"ksp_type": "preonly", "pc_type": "lu","pc_factor_mat_solver_type": "mumps"})

#solve adj problem in lamda (lamda is adjoint variable)
lamda = problem_adj.solve()

print("The objective is",J_old)

lamda_ = lamda.x.array
min_ = min(lamda_)
max_ = max(lamda_)
print("lamda min :",min_)
print("lamda max :",max_)

# Update dJde (sensitivity)
dJde = fem.Function(V)
grad_c = grad(c)
grad_lamda = grad(lamda)

dJde = (k0*gamma*((1-eps.x.array)**(eta-1))*c.x.array) + (D0*eta*(eps.x.array**(eta-1)))*(grad_c[0]*grad_lamda[0]+grad_c[1]*grad_lamda[1])+(k0*gamma*((1-eps.x.array)**(eta-1))*c.x.array*lamda.x.array)

print("dJde :",dJde)

The results are

The objective is (0.07877750212780599+0j)
lamda min : (-0.9155314560137311+0j)
lamda max : 0j
dJde : [Sum(FloatValue(0.07740840145414743), Sum(FloatValue(-0.08455854979804014), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
 Sum(FloatValue(0.07737913038121018), Sum(FloatValue(-0.08452332241286431), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
 Sum(FloatValue(0.07737327161484303), Sum(FloatValue(-0.08451627181655852), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
 ...
 Sum(FloatValue(-1.0), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),)))))))
 Sum(FloatValue(0.030059484951120978), Sum(FloatValue(-0.9689781605244571), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
 Sum(FloatValue(-1.0), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),)))))))]

I had also tried this method(You can see more explanations at this link Calculate sensitivity equation of reaction-diffusion systems with dolfinx - #20 by Pasuta) but the answer was also wrong and I think there is something wrong with this line.

dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))

# dJde is objactive function / eps
dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))
dJde = dolfinx.fem.petsc.assemble_vector(dJde_form)
dJde.assemble()

# dRde is the equation accounting for the transport of the species / eps
dRde_form = fem.form(adjoint(derivative(R,eps)))
dRde = dolfinx.fem.petsc.assemble_matrix(dRde_form)
dRde.assemble()

B = dJde.copy()
dRde.mult(lamda.vector, B)

# Update dJde
dJde += B
dJde.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

I will be very pleased. If someone can help solve this problem.

dokken · April 12, 2024, 11:45am

Here you are mixing numpy arrays (.x.array with ufl objects.

Pasuta:

had also tried this method(You can see more explanations at this link Calculate sensitivity equation of reaction-diffusion systems with dolfinx - #20 by Pasuta) but the answer was also wrong and I think there is something wrong with this line.

dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))

Please elaborate on what you imply by «wrong». This gives you the symbolic expression for dJ/d(eps), which you in turn can assemble into a vector

Pasuta · April 16, 2024, 3:09pm

As you can see the equation shown in the picture ,the first term equation is the same term as dJde in this line

dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))
dJde1 = dolfinx.fem.petsc.assemble_vector(dJde_form)
dJde1.assemble()
print("dJde1:",dJde1.array)

dJde1: [-4.22792750e-06+0.j -1.40883944e-06+0.j -4.22651812e-06+0.j ...
 -4.94829826e-05+0.j -4.87114826e-05+0.j -1.65374090e-05+0.j]

If I try substituting the values into the equation directly (in first term only), the variable is defined as follows and the maximum value of c is 1, so the maximum value of djde1 should be -1, But the maximum value printed is -1.65374090e-05 So I think that line of code is incorrect.

#Parameter
D0 = 0.1
k0 = -1.0
eta = 2.0
gamma = 2.0
eps = 0.5

And for this, I tried so many ways but still couldn’t find the correct way, can you correct it?? and I don’t know what is different with numpy array and ufl objects.

dJde = (k0gamma((1-eps.x.array)(eta-1))c.x.array) + (D0eta*(eps.x.array(eta-1)))*(grad_c[0]grad_lamda[0]+grad_c[1]grad_lamda[1])+(k0gamma((1-eps.x.array)**(eta-1))c.x.arraylamda.x.array)

dokken · June 23, 2024, 7:13pm

I don’t understand why the maximum value should be -1. Could you please elaborate on this?

dJde = -1 * 2 * (0.5)*c
with your current settings, so the first term should be dJde=-c?

If you are talking about:

then this is a very different derivative than the first one, as it is an integral over the domain, which is not the same as the values at the dofs.

Topic		Replies	Views
Calculate sensitivity equation of reaction-diffusion systems with dolfinx General mesh , assembly , dolfinx	22	629	April 9, 2024
Topology optimization loop cannot update the new design variable General optimization , dolfinx	11	292	March 14, 2024
How to set Neumann boundary conditions to solve adjoint variable variational formulation boundary-conditions , dolfinx	6	274	March 27, 2024
Error with complex numbers in dolfinx dolfinx	15	1415	August 19, 2022
Adjoint method using nonlinear residual and mixed elements dolfinx	1	194	March 26, 2024

Try to evaluate the sensitivity value but couldn't solve the correct value

Related topics