Hi,
I was wondering whether in dolfinx was possible to write the weak form below,
\int_{\Omega} \frac{\partial v}{\partial x_{i}} \, \frac{\partial u}{\partial x_{i}} \, dV - \, \int_{\Omega} v \, Q \, dV \, = 0
such as Q has already a prescribed size, say as a matrix of size [N,N].
Thank you.
Your question is unclear. For instance, what does N represent?
thanks for the reply. In this context, N would mean number of surfaces of part of my mesh. Is it clear now?
Where in your matrix would you like to insert it? Will it be a dense matrix?
How does the different N entries depend on eachother?
Finally, do you want this to work in parallel? If so, could you give a minimal reproducible example of what you have so far?
thank you for the replying my post. I did not understand the question âWhere in your matrix would you like to insert it?â. For the other questions:
. âWill it be a dense matrix?â: Yes;
. âHow does the different N entries depend on eachother?â: They should be constants, not being affected by the state variables of my PDE. They have to be organized in a proper fashion though, inside the matrix [Q]. But that would be my duty, I guess;
. âFinally, do you want this to work in parallel?â: Yes.
Below follows a mwe, that generates the plot we see in the figure. The idea here would be to, instead of having a single value for âheat_flux_valueâ variable, I would have a matrix of, in this example, [10x10] values. Iâm using dolfinx 0.6.0:
from mpi4py import MPI
from dolfinx.fem import (
Function, FunctionSpace, form, locate_dofs_topological, dirichletbc, Constant
)
from dolfinx.fem.petsc import (
apply_lifting, assemble_matrix, assemble_vector, create_vector, set_bc
)
from dolfinx.mesh import (
CellType, create_box, locate_entities_boundary, meshtags, locate_entities
)
from ufl import (
dx, grad, inner, TrialFunction, TestFunction, Measure
)
from petsc4py.PETSc import ScalarType
from dolfinx.io import XDMFFile
import numpy as np
from dolfinx import fem, cpp, log
from petsc4py import PETSc
k = 1.0
heat_flux_value = 1e4
lx, ly, lz = 2.0, 1.0, 1.0
nelx, nely, nelz = 20, 10, 10
number_element = nelx*nely*nelz
domain = create_box(comm=MPI.COMM_WORLD, points=[(0., 0., 0.), (lx, ly, lz)], n=[nelx, nely, nelz], cell_type=CellType.hexahedron)
H = FunctionSpace(domain, ("CG", 1))
phi = TrialFunction(H)
h = TestFunction(H)
temperature = Function(H)
D = FunctionSpace(domain, ("DG", 0))
tolerance = 1e-6
boundary_conditions_heat_transfer = []
def temperature_zero(x):
return x[0] > (lx - tolerance)
facets_zero_temperature = locate_entities_boundary(domain, domain.topology.dim-1, temperature_zero)
zero_temperature = Function(H)
zero_temperature.x.set(0)
bcu_front_temperature = dirichletbc(zero_temperature,\
locate_dofs_topological(H, domain.topology.dim-1,\
facets_zero_temperature))
boundary_conditions_heat_transfer = [bcu_front_temperature]
def heat_flux(x):
return x[0] < tolerance
facets_thermal_load = locate_entities_boundary(domain, domain.topology.dim-1, heat_flux)
facet_indices_thermal_load = [facets_thermal_load]
facet_markers_thermal_load = []
facet_markers_thermal_load.append(np.full_like(facet_indices_thermal_load[0], 1))
facet_indices_thermal_load = np.hstack(facet_indices_thermal_load).astype(np.int32)
facet_markers_thermal_load = np.hstack(facet_markers_thermal_load).astype(np.int32)
sorted_facets_thermal_load = np.argsort(facet_indices_thermal_load)
facet_tags_thermal_load = meshtags(domain, domain.topology.dim-1,\
facet_indices_thermal_load[sorted_facets_thermal_load],\
facet_markers_thermal_load[sorted_facets_thermal_load])
phi = TrialFunction(H)
h = TestFunction(H)
dx = Measure("dx", domain=domain, metadata={"quadrature_degree": 2})
dOmega = Measure("ds", domain=domain, subdomain_data=facet_tags_thermal_load, metadata=None)
bilinear = k*inner(grad(h),grad(phi))*dx
linear = heat_flux_value*h*dOmega(1)
bilinear_form = form(bilinear)
linear_form = form(linear)
temperature = Function(H)
bilinear_matrix = assemble_matrix(bilinear_form, boundary_conditions_heat_transfer)
bilinear_matrix.assemble()
linear_vector = create_vector(linear_form)
assemble_vector(linear_vector, linear_form)
apply_lifting(linear_vector, [bilinear_form], [boundary_conditions_heat_transfer])
linear_vector.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(linear_vector, boundary_conditions_heat_transfer)
options = PETSc.Options()
linear_solver = PETSc.KSP().create(MPI.COMM_WORLD)
linear_solver.setType("cg")
linear_solver.setTolerances(rtol=1e-6)
linear_solver.getPC().setType("gamg")
options["ksp_monitor"] = None
linear_solver.setFromOptions()
linear_solver.setInitialGuessNonzero(False)
linear_solver.setOperators(bilinear_matrix)
linear_solver.solve(linear_vector, temperature.vector)
temperature.x.scatter_forward()
with XDMFFile(domain.comm, "temperature.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(temperature)
Please format the code with 3x`
Ie
```python
# add code here
```
Unfortunately,
I do still not understand how an index in Q relates to a position in your matrix.
It is also unclear to mean how the latter matrix should be assembled, as V is a scalar function space.
Hi @dokken,
below follows a figure where I tried to show you how the heat fluxes matrix is organized in respect to the position of each of its value in a 10x10 mesh.
Thank you.
A quick and dirty way would be something like
linear = h*dOmega(1)
linear_form = form(linear)
linear_vector = create_vector(linear_form)
assemble_vector(linear_vector, linear_form)
Q_matrix_np = # your matrix in numpy format
linear_vector.array[:] = Q_matrix_np @ linear_vector.array
It wonât work in parallel, but at least should give you an idea about how to start.
Iâm still confused about the question overall. In particular, if Q really is dense, that this isnât like any other heat flux model Iâve ever seen, because having Q dense means that the temperature in a point of the domain affects the temperature everywhere else, and not only in the immediate vicinity.
Q is dense only at the face where heat flux is prescribed, like in the figure below:
My application demands to have a particular distribution of heat fluxes (known beforehand) over that face, thatâs the origin of my question. So, per your answer, am I able to embed a vector of prescribed values created through âcreate_vectorâ in the weak form?
You can give a try at the snippet I wrote above.
