Hi everybody,
I’m a FEniCS and python newbie, and I am trying to implement a small program
that solves Gauss’s law.
I already have something works reasonable well in a uniform mesh, and I am
trying to implement adaptive meshes.
I have several problems on this respect.
The first one is with the AdaptiveLinearSolve.
I am using the Energy as the functional to be optimised (i.e. E=0.5*rho*phi,
rho being the charge distribution, phi being the potential calculated with
Linear solver)
Please find the minimal portion of the script attached.
The problem is, the energy calculated by the assemble, and reported in the
solver.summary() are very different.
Also, rho being a normalised gaussian charge distribution, it should integrate
to 1 (which it does in an uniform mesh with converged lattice space), however,
after AdaptiveLinearSolve, the integral goes to 0.1 or some other strange
value.
I suspect the Function(V) does not use the optimised mesh but the original
coarse mesh. Can it be the case? What else can it be? How can it be corrected?
Best,
Baris
Minimal portion of the script (celldm,alat,MeshSize are predetermined
constants; epsilon,rho are predetermined functions; DirichletXYZ is a python
class constructed using “near” and on_boundary)
mesh =
BoxMesh(0.0,0.0,0.0,celldm1*alat,celldm2*alat,celldm3*alat,MeshSize,MeshSize,MeshSize)
V = FunctionSpace(mesh, "CG", 1)
u0 = Constant(0.0)
db =
DrichletXYZ(xmin=0.0,xmax=celldm1*alat,ymin=0.0,ymax=celldm2*alat,zmin=0.0,zmax=celldm3*alat)
bc = DirichletBC(V, u0, db)
bcs = [bc]
phi = TrialFunction(V)
psi = TestFunction(V)
a = (epsilon*inner(grad(phi), grad(psi)))*dx
L = 4*pi*rho*psi*dx
u = Function(V)
E = 0.5*u*rho*dx
problem = LinearVariationalProblem(a, L, u,bc)
solver = AdaptiveLinearVariationalSolver(problem, E)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"]
= "gmres"
solver.parameters["linear_variational_solver"]["linear_solver"] = "gmres"
solver.solve(MeshTol)
solver.summary()
# Report results
chg_ongrid=Function(V)
chg_ongrid.interpolate(rho)
chg = chg_ongrid*dx
chg_int = assemble(chg)
#E integral
E = u*rho*dx
E_int = assemble(E)
E_int = E_int*0.5
_______________________________________________
fenics-support mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics-support
