Regarding charge distribution being ~0.1: This was a due to a python mistake I made. (conversion to float was silently failing in a variable)
The problem with the Energy term being different continues. 28 Aug 2014 tarihinde 10:49 saatinde, obm <[email protected]> şunları yazdı: > 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
