Assertion in 0.6.0 which works fine in 0.5.0

marcelclaro · March 9, 2023, 5:28pm

I had a code working fine on version 0.5.0 which stopped to work in the new release 0.6.0, and I cannot find the problem. The subdomains are well defined.

The error:

Traceback (most recent call last):
  File "/home/marcel/KCMmodel/FourierIso.py", line 238, in <module>
    problem = LinearProblem(a, L, bcs=bc,  petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
  File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 566, in __init__
    self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options)
  File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 176, in form
    return _create_form(form)
  File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 171, in _create_form
    return _form(form)
  File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 138, in _form
    assert all([id(d) == id(data[0]) for d in data])
AssertionError

The code:

import numpy as np
import ufl
from dolfinx.fem.petsc import LinearProblem
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.io import XDMFFile
from petsc4py.PETSc import ScalarType
import gmsh
import sys
import argparse


#import system parameters
side = 40.0 #xy size
height = 100.0 #z size
dL = 10   #xy  general resolution
dH = 5   #z layer resolution
substrate_res = 20 #z substrate resolution

thickness = [0.06]  #Thickness of each layer
c=[2.49*1e-12,2.98*1e-12] #Specific heat c in J/(um)3K
k=[66*1e-6,33*1e-6] #isotropic conductivity in J/(um)2K
tbc = [5e-5]  #Thermal boundary conductance of each interface (same order)

pen_depth = 0.012

rpump = 5.0 #Pump beam radius
rprobe = 5.0 #Probe beam radius
beamoffset = 0.0 #Offset between paump and probe beam center

bumppar = -5 #Parameter to increase the resolution in the center
progressionpar = 1.2 #Parameter to increase the resolution close to the surface (1.1-1.3 probably works)


#calculate interface z-coordinates
interfaces = np.cumsum(np.array(thickness))
interfaces = np.append(interfaces,height)

#The frequency "omega" comes from the parameter
parser = argparse.ArgumentParser(description='Calculate phase shift.')
parser.add_argument('freq', type=float, nargs=1,
                    help='Frequency in Hz')
args = parser.parse_args()
omega = 2*np.pi*args.freq[0]


proc = MPI.COMM_WORLD.rank #For MPI control


# create the mesh
if proc == 0:

    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal",0)
    gmsh.model.add("Model")
    lc = 10 #point resolution - not really used
    gmsh.model.geo.addPoint(0, 0, 0, lc, 1) #layer 1 points
    gmsh.model.geo.addPoint(side, 0, 0, lc, 2)
    gmsh.model.geo.addPoint(side, side, 0, lc, 3)
    gmsh.model.geo.addPoint(0, side, 0, lc, 4)

    gmsh.model.geo.addLine(1, 2, 1) #border layer 1
    gmsh.model.geo.mesh.setTransfiniteCurve(1, dL, "Bump", bumppar)
    gmsh.model.geo.addLine(2, 3, 2)
    gmsh.model.geo.mesh.setTransfiniteCurve(2, dL, "Bump", bumppar)
    gmsh.model.geo.addLine(3, 4, 3)
    gmsh.model.geo.mesh.setTransfiniteCurve(3, dL, "Bump", bumppar)
    gmsh.model.geo.addLine(4, 1, 4)
    gmsh.model.geo.mesh.setTransfiniteCurve(4, dL, "Bump", bumppar)

    gmsh.model.geo.addCurveLoop([-1, -2, -3, -4], 1)
    gmsh.model.geo.addPlaneSurface([1], 1) #1st surface
    gmsh.model.geo.mesh.setTransfiniteSurface(1)

    gmsh.model.geo.synchronize()

    leftwall_list = []
    rightwall_list = []
    frontwall_list = []
    backwall_list = []
    interfaces_index = []

    for count, i_height in enumerate(interfaces):
        gmsh.model.geo.addPoint(0, 0, i_height, lc, 5+4*count)
        gmsh.model.geo.addPoint(side, 0, i_height, lc, 6+4*count)
        gmsh.model.geo.addPoint(side, side, i_height, lc, 7+4*count)
        gmsh.model.geo.addPoint(0, side, i_height, lc, 8+4*count)

        gmsh.model.geo.addLine(5+4*count, 6+4*count, 5+8*count)
        gmsh.model.geo.mesh.setTransfiniteCurve(5+8*count, dL, "Bump", bumppar)
        gmsh.model.geo.addLine(6+4*count, 7+4*count, 6+8*count)
        gmsh.model.geo.mesh.setTransfiniteCurve(6+8*count, dL, "Bump", bumppar)
        gmsh.model.geo.addLine(7+4*count, 8+4*count, 7+8*count)
        gmsh.model.geo.mesh.setTransfiniteCurve(7+8*count, dL, "Bump", bumppar)
        gmsh.model.geo.addLine(8+4*count, 5+4*count, 8+8*count)
        gmsh.model.geo.mesh.setTransfiniteCurve(8+8*count, dL, "Bump", bumppar)

        gmsh.model.geo.addCurveLoop([(5+8*count), (6+8*count), (7+8*count), (8+8*count)], 2+5*count) #top surface edge loops

        if count == len(interfaces)-1:
            gmsh.model.geo.addLine(1+4*count, 5+4*count, 9+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(9+8*count, substrate_res, "Progression", progressionpar)
            gmsh.model.geo.addLine(2+4*count, 6+4*count, 10+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(10+8*count, substrate_res, "Progression", progressionpar)
            gmsh.model.geo.addLine(3+4*count, 7+4*count, 11+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(11+8*count, substrate_res, "Progression", progressionpar)
            gmsh.model.geo.addLine(4+4*count, 8+4*count, 12+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(12+8*count, substrate_res, "Progression", progressionpar)
        else:
            gmsh.model.geo.addLine(1+4*count, 5+4*count, 9+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(9+8*count, dH)
            gmsh.model.geo.addLine(2+4*count, 6+4*count, 10+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(10+8*count, dH)
            gmsh.model.geo.addLine(3+4*count, 7+4*count, 11+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(11+8*count, dH)
            gmsh.model.geo.addLine(4+4*count, 8+4*count, 12+8*count)
            gmsh.model.geo.mesh.setTransfiniteCurve(12+8*count, dH)

        if count == 0:
            gmsh.model.geo.addCurveLoop([1, 10, -5, -9], 3)
            gmsh.model.geo.addCurveLoop([2, 11, -6, -10], 4)
            gmsh.model.geo.addCurveLoop([3, 12, -7, -11], 5)
            gmsh.model.geo.addCurveLoop([4, 9, -8, -12], 6)
        else:
            gmsh.model.geo.addCurveLoop([(8*count-3), (10+8*count), -(5+8*count), -(9+8*count)], 3+5*count)
            gmsh.model.geo.addCurveLoop([(8*count-2), (11+8*count), -(6+8*count), -(10+8*count)], 4+5*count)
            gmsh.model.geo.addCurveLoop([(8*count-1), (12+8*count), -(7+8*count), -(11+8*count)], 5+5*count)
            gmsh.model.geo.addCurveLoop([(8*count), (9+8*count), -(8+8*count), -(12+8*count)], 6+5*count)

        gmsh.model.geo.addPlaneSurface([2+5*count], 6+5*count) #top surface is the last
        gmsh.model.geo.mesh.setTransfiniteSurface(6+5*count)
        
        gmsh.model.geo.addPlaneSurface([3+5*count], 3+5*count)
        gmsh.model.geo.mesh.setTransfiniteSurface(3+5*count)
        gmsh.model.geo.addPlaneSurface([4+5*count], 4+5*count)
        gmsh.model.geo.mesh.setTransfiniteSurface(4+5*count)
        gmsh.model.geo.addPlaneSurface([5+5*count], 5+5*count)
        gmsh.model.geo.mesh.setTransfiniteSurface(5+5*count)
        gmsh.model.geo.addPlaneSurface([6+5*count], 2+5*count)
        gmsh.model.geo.mesh.setTransfiniteSurface(2+5*count)

        
        lastsurface = 6+5*count #keep track of last surface
        interfaces_index.append(6+5*count)
        frontwall_list.append(3+5*count)
        rightwall_list.append(4+5*count)
        backwall_list.append(5+5*count)
        leftwall_list.append(2+5*count)
        
        #create layer volume
        gmsh.model.geo.addSurfaceLoop([1+5*count, 2+5*count, 3+5*count, 4+5*count, 5+5*count, 6+5*count])
        gmsh.model.geo.addVolume([1+count], 1+count)
        gmsh.model.geo.mesh.setTransfiniteVolume(1+count)
        gmsh.model.geo.synchronize()
        gmsh.model.addPhysicalGroup(3, [1+count], 1+count)

    interfaces_index.pop(-1)
    gmsh.model.addPhysicalGroup(2, [lastsurface], 2)
    gmsh.model.addPhysicalGroup(2, [1], 1)
    gmsh.model.addPhysicalGroup(2, frontwall_list, 3)
    gmsh.model.addPhysicalGroup(2, rightwall_list, 4)
    gmsh.model.addPhysicalGroup(2, backwall_list, 5)
    gmsh.model.addPhysicalGroup(2, leftwall_list, 6)
    for n, i in enumerate(interfaces_index):
        gmsh.model.addPhysicalGroup(2, [i], 7+n)
    gmsh.model.geo.synchronize()
 
    #gmsh.option.setNumber("Mesh.RecombineAll", 1)
    gmsh.model.mesh.generate(3)

    gmsh.option.setNumber("Mesh.SaveAll", 0)
    gmsh.model.geo.synchronize()
    gmsh.write("tempmesh.msh")

    #if '-nopopup' not in sys.argv:
        #gmsh.fltk.run()
    gmsh.finalize()

#read the mesh file
from dolfinx.io.gmshio import read_from_msh
model_rank = 0
domain, cell_tags, facet_tags = read_from_msh("tempmesh.msh", MPI.COMM_WORLD, model_rank, gdim=3)

#Create function space
V = fem.FunctionSpace(domain, ("DG", 1)) #Scalar field works for Fourier
x = ufl.SpatialCoordinate(domain) #spatial coordinates

#define the heat flux at the surface and the volumetric heating
f =  2e-6*((ufl.exp(-x[2]/pen_depth))) * x[0] * (x[0]-side) * x[1] * (x[1]-side) * (2 / (np.pi * rpump * rpump) )  * ufl.exp( - 2 * ((x[0]-side/2.0)**2 + (x[1]-side/2.0)**2) / (rpump * rpump ))


#Variatinol problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

#Find boundary points
def boundary_one(x):
    return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0],side))
def boundary_two(x):
    return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1],side))
def boundary_three(x):
    return np.isclose(x[2], height)

#define dirichlet boundary conditions
u_D = fem.Function(V, dtype=np.complex128)
u_D.interpolate(lambda x: 0.0 + 0.0j + x[0]*0 + x[1]*0 + x[2]*0)
bc = []
dofs_D = fem.locate_dofs_geometrical(V, boundary_one)
bc.append(fem.dirichletbc(u_D, dofs_D))
dofs_D = fem.locate_dofs_geometrical(V, boundary_two)
bc.append(fem.dirichletbc(u_D, dofs_D))
dofs_D = fem.locate_dofs_geometrical(V, boundary_three)
bc.append(fem.dirichletbc(u_D, dofs_D))

#Create functions for k and c in the space
Q = fem.FunctionSpace(domain, ("DG", 0))
kappa = fem.Function(Q)
heatcap = fem.Function(Q)
for layer in range(0,len(interfaces)):
    foundcells = cell_tags.find(layer+1)
    kappa.x.array[foundcells] = np.full_like(foundcells, k[layer], dtype=ScalarType)
    heatcap.x.array[foundcells] = np.full_like(foundcells, c[layer], dtype=ScalarType)


#tag surface for selective integration
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_tags)

mu = 1 #this is the internal penalty parameter, 1e-2 works well so far
n = ufl.FacetNormal(domain) #normal of each surface

#Fourier problem using SIPM, as  10.1016/j.cam.2006.08.029
a  = ( kappa*ufl.inner(ufl.grad(u), ufl.grad(v))  + 1j* omega * heatcap * ufl.inner(u,v) ) * ufl.dx  \
    - ufl.conj(ufl.inner(ufl.avg(kappa*ufl.grad(v)),ufl.jump(u,n))) * ufl.dS \
    - ufl.inner(ufl.avg(kappa*ufl.grad(u)),ufl.jump(v,n))*ufl.dS \
    + mu*ufl.inner(ufl.jump(u),ufl.jump(v)) * ufl.dS
#For each interface is necessary to remove general interface conditions and add the interface boundary conductivity
for i, h in enumerate(tbc):
    a -= mu*ufl.inner(ufl.jump(u),ufl.jump(v))*dS(7+i)
    a += ufl.inner(ufl.jump(u,n)*h,ufl.jump(v,n))*dS(7+i)
    a += ufl.conj(ufl.inner(ufl.avg(kappa*ufl.grad(v)),ufl.jump(u,n)))*dS(7+i)
    a += ufl.inner(ufl.avg(kappa*ufl.grad(u)),ufl.jump(v,n))*dS(7+i)
L =   + ufl.inner(f, v) * ufl.dx 


problem = LinearProblem(a, L, bcs=bc,  petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

#Export results to a file
import dolfinx.io
uh.name = "Temperature"
kappa.name = "k"
heatcap.name = "c"
with dolfinx.io.XDMFFile(domain.comm, "outputfourier.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(uh)
    xdmf.write_function(kappa)
    xdmf.write_function(heatcap)

dokken · March 9, 2023, 8:34pm

Your error is not reproducible, as you have not supplied

and a whole lot of parameters.

The check that is failing in your code is:

github.com

FEniCS/dolfinx/blob/v0.6.0/python/dolfinx/fem/forms.py#L136-L139


      
          # Check that the subdomain data for each integral of this type is
          # the same
          assert all([id(d) == id(data[0]) for d in data])
          subdomains[integral_type] = data[0]

which checks that you have supplied the same subdomain-data (meshtag/facettag) for every integral of the same kind.

As I cannot reproduce your error (as stated above), I cannot go into more detail on what could have gone wrong in your problem.

marcelclaro · March 10, 2023, 8:10am

Thank you dokken.
I edited the code to include those parameters.

I think the problem was mix ufl.dS with the defined measure. It seems do work now.

#tag surface for selective integration
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_tags)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_tags)

mu = 1 #this is the internal penalty parameter, 1e-2 works well so far
n = ufl.FacetNormal(domain) #normal of each surface

#Fourier problem using SIPM, as  10.1016/j.cam.2006.08.029
a  = ( kappa*ufl.inner(ufl.grad(u), ufl.grad(v))  + 1j* omega * heatcap * ufl.inner(u,v) ) * dx  \
    - ufl.conj(ufl.inner(ufl.avg(kappa*ufl.grad(v)),ufl.jump(u,n))) * dS \
    - ufl.inner(ufl.avg(kappa*ufl.grad(u)),ufl.jump(v,n))*dS \
    + mu*ufl.inner(ufl.jump(u),ufl.jump(v)) * dS
#For each interface is necessary to remove general interface conditions and add the interface boundary conductivity
for i, h in enumerate(tbc):
    a -= mu*ufl.inner(ufl.jump(u),ufl.jump(v))*dS(7+1)
    a += ufl.inner(ufl.jump(u,n)*h,ufl.jump(v,n))*dS(7+1)
    a += ufl.conj(ufl.inner(ufl.avg(kappa*ufl.grad(v)),ufl.jump(u,n)))*dS(7+1)
    a += ufl.inner(ufl.avg(kappa*ufl.grad(u)),ufl.jump(v,n))*dS(7+1)
L =   + ufl.inner(f, v) * dx

Thank you

Topic		Replies	Views
AttributeError when using dolfinx.fem.petsc.LinearProblem Errors petsc	1	57	January 16, 2025
Assertion error, mesh is None General	1	291	August 5, 2023
Assertion error on Poisson problem with complex option General	9	483	March 25, 2024
AttributeError: module 'dolfinx.cpp.la.petsc' has no attribute 'create_vector_wrap' dolfinx	1	232	January 12, 2024
Fem.petsc not found in complex-mode dolfinx petsc	6	370	June 15, 2024

Assertion in 0.6.0 which works fine in 0.5.0

Related topics