Computing fluxes accurately

mathematics
21

Hi,
this works exactly the same, independently of mechanics. If you take the action of the thermal bilinear form on a field with unit value on a given boundary, you will get the total flux on this boundary in a consistent manner.

1 Like
22

Thanks a lot Jeremy, I am very excited, I will post my (failed) attempt when I have the time. So far my whole personal project solely depends on evaluating these fluxes as accurately as possible, I was thinking about ways to bypass the problem but I did not find any remedy. I now have big hopes!

23

My current attempt looks like so:


import numpy as np
from dolfinx import log, fem
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector, action)
import gmsh
from dolfinx.io import gmshio, VTXWriter


mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/mesh2.msh", MPI.COMM_WORLD, gdim=2)

def thermal_problem(ΔT=0):
    current_curves=(1,2)
    voltages_curves=(1,2)  
    inj_current_curve, out_current_curve = current_curves
    reading_voltage_curve_0, reading_voltage_curve_1 = voltages_curves

    # Define FE function space
    deg = 1
    el = FiniteElement("CG", mesh.ufl_cell(), deg)
    V = FunctionSpace(mesh, el)    

    u = TestFunction(V)
    temp = Function(V)
  
    κ = 1.8
    T_cold = 300
    # Define the boundary conditions
    left_facets = facet_markers.find(inj_current_curve)
    right_facets = facet_markers.find(out_current_curve)
  
    left_dofs_temp = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
    right_dofs_temp = locate_dofs_topological(V, mesh.topology.dim-1, right_facets)
    
    bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, V)
    bc_temp_right = dirichletbc(ScalarType(T_cold + ΔT), right_dofs_temp, V)
    bcs = [bc_temp_left, bc_temp_right]
    
    x = SpatialCoordinate(mesh)
    dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
    ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)


    # Weak form.
    weak_form = dot(-κ * grad(temp), grad(u)) * dx

    print(f''' ------- Pre-processing --------
    Length of the side where heat enters: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
    Length of the side where heat leaves the material: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
    ''')

    # Solve the PDE.
    problem = NonlinearProblem(weak_form, temp, bcs=bcs)
    solver = NewtonSolver(MPI.COMM_WORLD, problem)

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}pc_type"] = "lu"
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"

    ksp.setFromOptions()

    log.set_log_level(log.LogLevel.WARNING)
    n, converged = solver.solve(temp)
    assert (converged)
    print(f'''------- Processing --------
    Number of interations: {n:d}''')

    avg_temp = assemble_scalar(form(temp * dx(domain=mesh))) / assemble_scalar(form(1 * dx(domain=mesh)))
    print(f'avg_temp: {avg_temp}')

    # Compute fluxes on boundaries
    n = FacetNormal(mesh)
    down_heat_flux = form(dot(-κ * grad(temp), n)*ds(inj_current_curve))
    Q_out = assemble_scalar(down_heat_flux)

    top_heat_flux = form(dot(-κ * grad(temp), n)*ds(out_current_curve))
    Q_in = assemble_scalar(top_heat_flux)

    print(f'''------- Post processing --------
    Q_in: {Q_in}
    Q_out: {Q_out}''')    


    Q = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
    q = Function(Q)
    flux_calculator = Expression(-κ * grad(temp), Q.element.interpolation_points())
    q.interpolate(flux_calculator)
    
    # Try an alternate way to compute the fluxes.
    virtual_temp = Function(V)
    virtual_heat_flux_form = form(action(weak_form, virtual_temp))
    virtual_heat_flux = assemble_scalar(virtual_heat_flux_form)
    print(f'Virtual heat flux computation: {virtual_heat_flux}')
    return Q_in, Q_out
    
t_b = thermal_problem(ΔT=10)
print(t_b)

Producing the output:

 ------- Pre-processing --------
    Length of the side where heat enters: 0.10000000000000002
    Length of the side where heat leaves the material: 0.09999999999999998
    
------- Processing --------
    Number of interations: 1
avg_temp: 305.4445644560497
------- Post processing --------
    Q_in: -7.434551004310087
    Q_out: 7.494950253817867
Virtual heat flux computation: 0.0
(-7.434551004310087, 7.494950253817867)

In other words, it computes as 0 heat flux instead of something like 7.49 or so.
My mesh is similar to the file.msh I posted in post #10 just above.
Any help is appreciated.