Error using spatially varying anisotropic diffusion coefficient

[bragostin]

Hello,

I am setting-up a 3D heat diffusion problem using spatially a varying anisotropic diffusion coefficient in the form:

    xc, yc, zc = baseMesh.cellCenters
    tyzc = (zc<=(yc-A)) & (zc>=(yc-B)) & (zc>=(yc+C)) & (zc<=(yc+D))
    t_tim = (xc < 0.) & (xc > -etim) & tyzc
    t_mb = (xc < -etim) & (xc > -(etim+emb)) & tyzc
    t_sd = (xc < -(etim+emb)) & (xc > -(etim+emb+esd)) & tyzc
    t_ms = (xc < -(etim+emb+esd)) & (xc > -(etim+emb+esd+ems)) & tyzc
    t_mt = (xc < -(etim+emb+esd+ems)) & (xc > -(etim+emb+esd+ems+emt)) & tyzc
    t_ins = (xc < 0.) & (xc > -(etim+emb+esd+ems+emt)) & N.logical_not(tyzc)
    k_x = k_x_tim*t_tim + kmb*t_mb + ksd*t_sd + kms*t_ms + kmt*t_mt + 1e-6*t_ins
    k_y = k_yz_tim*t_tim + kmb*t_mb + ksd*t_sd + kms*t_ms + kmt*t_mt + 1e-6*t_ins
    k_z = k_yz_tim*t_tim + kmb*t_mb + ksd*t_sd + kms*t_ms + kmt*t_mt + 1e-6*t_ins
    Dk = CellVariable(mesh=baseMesh, value=[[k_x,0.,0.],[0.,k_y,0.],[0.,0.,k_z]])
    equation = ImplicitDiffusionTerm(Dk) + ImplicitSourceTerm(boundarySource.divergence) - (FrontFaces * frontValue).divergence

It represents a anisotropic diffusion coefficient varying spatially through different layers of materials. The anisotropy is present only in one of the layers (k_x_tim, k_yz_tim)

I get the following error message when solving equation with equation.solve(T). Everything works fine when I am using an isotropic diffusion coefficient. Any idea?

 File "C:\WinPython-64bit-3.4.4.4Qt5\python-3.4.4.amd64\lib\site-packages\scipy\sparse\sputils.py", line 51, in upcast
   raise TypeError('no supported conversion for types: %r' % (args,))
TypeError: no supported conversion for types: (dtype('O'),)

[bragostin]

Some update, I found something that works:

Dk = FaceVariable(mesh=baseMesh, rank=3, value=((1.,0.,0.),(0.,1.,0.),(0.,0.,1.)))
t_tim = (xc < 0.) & (xc > -etim) & tyzc
Dk[:, 0, :, t_tim] = k_x * Dk[:, 0, :, t_tim]
Dk[:, 1, :, t_tim] = k_yz * Dk[:, 1, :, t_tim]
Dk[:, 2, :, t_tim] = k_yz * Dk[:, 2, :, t_tim]

etc... for the different layers, and

equation = ImplicitDiffusionTerm(Dk[0]) + ImplicitSourceTerm(boundarySource.divergence) - (FrontFaces * frontValue).divergence
equation.solve(T)

No more error message. However, I tried to validate with the isotropic case by comparing

Dk = FaceVariable(mesh=baseMesh, value=1.)
Dk.setValue(k_x, where=(x < 0) & (x > -etim) & tyz)
equation = ImplicitDiffusionTerm(coeff=Dk) + ImplicitSourceTerm(boundarySource.divergence) - (FrontFaces * frontValue).divergence

which is the isotropic case with

Dk = FaceVariable(mesh=baseMesh, rank=3, value=((1.,0.,0.),(0.,1.,0.),(0.,0.,1.)))
t_tim = (xc < 0.) & (xc > -etim) & tyzc
Dk[:, :, :, t_tim] = k_x
equation = ImplicitDiffusionTerm(Dk[0]) + ImplicitSourceTerm(boundarySource.divergence) - (FrontFaces * frontValue).divergence

which is an anisotropic setting with equal diffusion coefficients in all directions. It gives me quite different results, like 10°C difference. Should not those two cases yield the same result?

[bragostin]

I think I found the solution:

t_tim = (x < 0.) & (x > -etim) & tyz
Dk = FaceVariable(mesh=baseMesh, rank=3, value=((1.,0.,0.),(0.,1.,0.),(0.,0.,1.)))
Dk[:, 0, :, t_tim] = k_x * Dk[:, 0, :, t_tim]
Dk[:, 1, :, t_tim] = k_x * Dk[:, 1, :, t_tim]
Dk[:, 2, :, t_tim] = k_x * Dk[:, 2, :, t_tim]

gives the exact same results as the reference isotropic case.
Then the anisotropic case is

Dk = FaceVariable(mesh=baseMesh, rank=3, value=((1.,0.,0.),(0.,1.,0.),(0.,0.,1.)))
Dk[:, 0, :, t_tim] = k_x * Dk[:, 0, :, t_tim]
Dk[:, 1, :, t_tim] = k_yz * Dk[:, 1, :, t_tim]
Dk[:, 2, :, t_tim] = k_yz * Dk[:, 2, :, t_tim]

This seems to give plausible results. Is that the correct way to do it?