Thank you for the message.
I need this for my bachelor thesis and am a bit pressed for time.
This formula is to be solved:
\begin{equation}
(\nabla \phi) . (\nabla \rho) - \frac{\rho ^2}{\epsilon_0} = 0
\end{equation}
(How can I insert latex code here?).
Here is the MWE:
#Mesh Erzeugung:(mit GMSH)
from mpi4py import MPI
from dolfinx import mesh
import warnings
warnings.filterwarnings("ignore")
import gmsh
gmsh.initialize()
from dolfinx import fem
from dolfinx import nls
from dolfinx import log
from dolfinx import io
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, create_vector, set_bc
from dolfinx.graph import create_adjacencylist
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (XDMFFile, cell_perm_gmsh, distribute_entity_data, extract_gmsh_geometry,
extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh)
from dolfinx.mesh import create_mesh, meshtags_from_entities
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import numpy as np
import matplotlib.pyplot as plt
import tqdm.notebook
from mpi4py import MPI
from petsc4py import PETSc
import scipy.constants as c
#Domain-Definieren:
L = 1.
H = 1.
c_x = 0.50
c_y = 0.50
r = 0.01
gdim = 2
rank = MPI.COMM_WORLD.rank
if rank == 0:
rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
if rank == 0:
Mesh_domain = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()
#Verbinden aller Mesh-Punkte
Mesh_domain_marker = 1
if rank == 0:
volumes = gmsh.model.getEntities(dim=gdim)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], Mesh_domain_marker)
gmsh.model.setPhysicalName(volumes[0][0], Mesh_domain_marker, "Mesh_domain")
inlet_marker, outlet_marker, wall_marker, obstacle_marker, wall_top = 2, 3, 4, 5, 6
inflow, outflow, walls, obstacle, top = [], [], [], [], []
if rank == 0: #Walls sind Ground, Obstical ist Leiter
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H/2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H/2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L/2, 0, 0]):
walls.append(boundary[1])
elif np.allclose(center_of_mass, [L/2, H, 0]):
top.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
gmsh.model.addPhysicalGroup(1, top, wall_top)
gmsh.model.setPhysicalName(1, wall_top, "top")
if rank == 0:
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 8)
gmsh.option.setNumber("Mesh.RecombineAll", 2)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
if MPI.COMM_WORLD.rank == 0:
# Mesh Geometrie
x = extract_gmsh_geometry(gmsh.model)
# Topologie für jedes Element vom Mesh
topologies = extract_gmsh_topology_and_markers(gmsh.model)
num_cell_types = len(topologies.keys())
cell_information = {}
cell_dimensions = np.zeros(num_cell_types, dtype=np.int32)
for i, element in enumerate(topologies.keys()):
properties = gmsh.model.mesh.getElementProperties(element)
name, dim, order, num_nodes, local_coords, _ = properties
cell_information[i] = {"id": element, "dim": dim, "num_nodes": num_nodes}
cell_dimensions[i] = dim
# Sort elements by ascending dimension
perm_sort = np.argsort(cell_dimensions)
# Broadcast cell type data and geometric dimension
cell_id = cell_information[perm_sort[-1]]["id"]
tdim = cell_information[perm_sort[-1]]["dim"]
num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
cell_id, num_nodes = MPI.COMM_WORLD.bcast([cell_id, num_nodes], root=0)
if tdim - 1 in cell_dimensions:
num_facet_nodes = MPI.COMM_WORLD.bcast( cell_information[perm_sort[-2]]["num_nodes"], root=0)
gmsh_facet_id = cell_information[perm_sort[-2]]["id"]
marked_facets = np.asarray(topologies[gmsh_facet_id]["topology"], dtype=np.int64)
facet_values = np.asarray(topologies[gmsh_facet_id]["cell_data"], dtype=np.int32)
cells = np.asarray(topologies[cell_id]["topology"], dtype=np.int64)
cell_values = np.asarray(topologies[cell_id]["cell_data"], dtype=np.int32)
else:
cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
cells, x = np.empty([0, num_nodes], np.int64), np.empty([0, gdim])
cell_values = np.empty((0,), dtype=np.int32)
num_facet_nodes = MPI.COMM_WORLD.bcast(None, root=0)
marked_facets = np.empty((0, num_facet_nodes), dtype=np.int64)
facet_values = np.empty((0,), dtype=np.int32)
# Create distributed mesh
ufl_domain = ufl_mesh_from_gmsh(cell_id, gdim)
gmsh_cell_perm = cell_perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]
mesh = create_mesh(MPI.COMM_WORLD, cells, x[:, :gdim], ufl_domain)
tdim = mesh.topology.dim
fdim = tdim - 1
local_entities, local_values = distribute_entity_data(mesh, fdim, marked_facets, facet_values)
mesh.topology.create_connectivity(fdim, tdim)
adj = create_adjacencylist(local_entities)
# Create DOLFINx MeshTags
ft = meshtags_from_entities(mesh, fdim, adj, np.int32(local_values))
ft.name = "Facet tags"
#Funktionspace und Vektorfunktionspace defenieren
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, s_cg1)
V_E= fem.VectorFunctionSpace(mesh, ("DG", 0))
fdim = mesh.topology.dim - 1
#Konstanten, Funktionen und Randbedingungen definieren:
phi_0 = 80000 #V
rho0_0 = 1E-6 #C/m^3
T_0 = 20 #°C
p_0 = 1013 #hPa
h = c_y
#Kontroll-Variablen
LoopBreakPoint = 100
It = 0
CheckerRho = False
tol = 1e-1
num = 0
#Initialisierung der RB für rho
rho_aend = rho0_0
rho_0 = Constant(mesh, PETSc.ScalarType(rho0_0))
u1_0 = Constant(mesh, PETSc.ScalarType(phi_0))
E_0 = 33.7E+5 #V/m
K = 0.024 #sqrt(m)
delta = ((273.15 + 20)/(273.15 + T_0))*(p_0/1013)
E_c = E_0*delta*(1+(K/np.sqrt(delta*r))) # V/m
print("E_c: ", E_c)
#Boundary conditions definieren
# Walls
u1_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
u11_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
wall_facets = ft.indices[ft.values == wall_marker]
bcu1_walls = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(V, fdim, wall_facets), V)
# Obstacle
obstacle_facets = ft.indices[ft.values == obstacle_marker]
bcu1_obstacle = dirichletbc(PETSc.ScalarType(14e4), locate_dofs_topological(V, fdim, obstacle_facets), V)
#Boundary conditions for equations:
bcu1 = [bcu1_obstacle, bcu1_walls] # Spannung am Leiter und Ränder
#Laplace
u1 = ufl.TrialFunction(V)
v1 = ufl.TestFunction(V)
k = Constant(mesh, PETSc.ScalarType(0))
f = k/c.epsilon_0
epsilon = c.epsilon_0
a = ufl.inner(ufl.grad(u1),ufl.grad(v1)) * ufl.dx
L = ufl.inner(f , v1) * ufl.dx
u1 = fem.Function(V)
problem1 = fem.petsc.LinearProblem(a, L, bcu1, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi = problem1.solve()
#Bestimmung der Randbedingung für rho
while not CheckerRho:
print('solution for density is analysed')
It = It+1
print(It, "iteration of the PDE-system")
if It > LoopBreakPoint :
print('more then ', LoopBreakPoint, ' iterations')
break
phi_old = (phi.x.array.real)
E= - ufl.grad(phi)
E_expr = fem.Expression(E, V_E.element.interpolation_points)
EE = fem.Function(V_E)
EE.interpolate(E_expr)
rho = fem.Function(V)
v2 = ufl.TestFunction(V)
F = - ufl.inner((rho**2)/c.epsilon_0,v2)*dx + ufl.inner(ufl.dot(ufl.grad(phi),ufl.grad(rho)),v2)*dx # Kontinutätsgleichung
rho_BC = dirichletbc(PETSc.ScalarType(rho_aend), locate_dofs_topological(V, fdim, obstacle_facets), V)
problem2 = fem.petsc.NonlinearProblem(F, rho, bcs=[rho_BC])
#nonlinearproblem solver
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem2)
solver.max_it = 100
solver.atol = 1e-5
solver.rtol = 1e-12
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(rho)
assert(converged)
print(f"Number of interations for Newton-solver: {n:d}")
#Poisson
u1 = ufl.TrialFunction(V)
v1 = ufl.TestFunction(V)
aP = ufl.inner(grad(u1), grad(v1))*dx
LP = ufl.inner(-rho/c.epsilon_0,v1)*dx
u1 = fem.Function(V)
problem3 = fem.petsc.LinearProblem(aP, LP, bcu1, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi = problem3.solve()
#new Checker for Rho
phi_current= (phi.x.array.real)
print('finished itteration ', It)
print ('Iterations for potential have ended. Finished in ', It ,' itterations')
if np.abs((phi_current - phi_old)).all() < tol:
print("1 max: ",(np.max(np.abs(EE.x.array))))
E= -ufl.grad(phi)
E_expr = fem.Expression(E, V_E.element.interpolation_points)
EE = fem.Function(V_E)
EE.interpolate(E_expr)
print("2 max: ",(np.max(np.abs(EE.x.array))))
print(np.abs((np.max(np.abs(EE.x.array))-E_c))/E_c, " < ", tol)
print("rho_aend: ", rho_aend)
if np.abs((np.max(np.abs(EE.x.array))-E_c)/E_c)< tol:
CheckerRho = True
print("!!! finish !!!")
else:
CheckerRho = False
rho_aend += rho_aend*0.6*(((np.max(np.abs(EE.x.array))-E_c))/((np.max(np.abs(EE.x.array))+E_c)))
The problem, as I see it, is only in the while loop.
