Heat equation with pure Neumann boundary conditions

[Ruslan]

Hello there,

I’m trying to solve Heat equation in 2D with pure Neuman boundary conditions. And I can’t get normal solution. The problem doesn’t seem so complicated, but there is no proper example/demo for it. At least I didn’t find it.
I used this as start point:
https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/neumann-poisson/python/documentation.html
But got a lot of issues with expressions shapes, which I can’t handle.

I suppose it wiil be helpful for more than just me, if some one could provide an example of what the fenics code should look like for a given problem.

Thank you!

[dokken]

You are looking at an outdated version of FEniCS (1.4.0). The demo for the latest version can be found here.

[Ruslan]

Thank you for response!
But I’ve already fixed all the differences between the new and the old demo on my own. And Poisson with pure Neumann BC runs well. Problems arise in dynamics (Heat) in mixed psace.

[dokken]

Without you posting the code that you have, and what error message you Get, there isn’t much that can be done from my side.

[Ruslan]

I solved the problem using subdomains, and not using mixed space.
Some custom ploting method is also provided.
Here is the code:

import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from IPython.display import clear_output
import numpy as np
import itertools
from dolfin import *

# Create classes for defining parts of the boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

nx = 64
ny = 64

# Define mesh
mesh = UnitSquareMesh(nx, ny)

t = 0.0 # initial time 
Time = 1.0         # final time
num_steps = 20     # number of time steps
dt = Time / num_steps # time step size

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define input data
a0 = Constant(1.0)
a1 = Constant(0.1)
g = Constant("0.0001")
f = Constant("4.")

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)

T_d = Expression('sin(pi*x[1])*exp( -pi*pi*t )*pow(sin(pi*x[0]) , 3)', degree=2, t=0)
T_old = interpolate(T_d, V)

# Define new measures associated with the exterior boundaries
ds = ds(subdomain_data = boundaries)

# Define variational form
F = (u*v*dx - T_old*v*dx
    + dt*inner(a1*grad(u), grad(v))*dx
    - dt*f*v*dx - dt*g*v*ds(1) - dt*g*v*ds(2) - dt*g*v*ds(3) - dt*g*v*ds(4)
    )

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)
u = Function(V)

# Solve problem
for n in range(num_steps):

    # Update current time
    t += dt   
    # Compute solution
    solve(a == L, u)
    


    #3D Plot of the Solution  
    %matplotlib inline
    x = np.linspace(0.0, 1.0, nx + 1)
    y = np.linspace(0.0, 1.0, ny + 1)
    
    X, Y = np.meshgrid(x, y) 
    
    vertex_values_u = u.compute_vertex_values()
    Z_vert = vertex_values_u.reshape((nx + 1,ny + 1 ))
        
    fig = plt.figure()
    ax = plt.axes(projection='3d')

    ax.plot_surface(X, Y, Z_vert, rstride=1, cstride=1,
                cmap='viridis', edgecolor='none')
    #ax.set_title('t = ' +str(round(t, 2)))
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('T')
    #ax.set_zlim3d(-1.0 , 1.0)
    
    plt.show()
    
    print('t = %.2f' % (t)) 
     
    clear_output(wait=True)
    
    T_old.assign(u)

Some criticism is welcome

[dokken]

As g is the same on all boundaries, you do not need subdomains and could simply use the variational formulation:
F = (u*v*dx - T_old*v*dx + dt*inner(a1*grad(u), grad(v))*dx - dt*f*v*dx - dt*g*v*ds)
Otherwise I have no specific comments as what you are doing is quite straight forward.

On the visualization side, you could use Paraview to visualize the data (with eiter XDMF or VTK/PVD files).

[Ruslan]

Ok, I see.
Thank you!

[Mahdi_king]

would you mind please sharing the problem please ?

[Ruslan]

It’s posted above already!

[Mahdi_king]

Thank you for your reply ,
The demo can not run in my system.
would you mind writing the entire engineering problem here?

suppose : a engine cylinder heat transfer including following condition…or a heat transfer in pipe is like beyond…

I need the complete problem include the details ,

warm regards

[dokken]

If you read the previous posts in this example, you can find this link to a fully documented example on the Poisson equation. You can find several other documented demos at https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/

[Ruslan]

You can get basic knowledge on the entire process by looking at Fenics tutorial https://fenicsproject.org/tutorial/. Though some things from there are depricated in current Fenics version, you can get info consistently. Mr. Dokkens second link provides acces to actualized demos from tutorial. Also as I mentioned before, I used subdomains. So this is the point, you might looking at.

[SHilbard]

Hi Ruslan
I have been looking at the same problem and used the same example as you as my starting point and I wonder why you did not include the constant that is in the example in your solution.
It is clearly easier to solve the problem with just one FunctionSpace instead of two but how do you avoid getting the error that, the constant is intended to account for ?

[Ruslan]

Hi,

Which constant do you mean? And why does it produce an error?

[SHilbard]

In the example they mention that we get an error because we only have Neumann conditions and therefor we can only know the flux in and out of the domain but not the value in the domain to begin with. And to account for this error we add a constant c.
And then when we write the variational from we get that we are looking for u and c with testfunctions v and d.

Can we just ignore that problem when we do time-dependent calculations because we define a starting value ?

[Ruslan]

I’m not an expert, but I would say yes.
See for instance:

Steady, nonlinear heat equation with pure neumann BCs

I believe there are two mistakes in the formulation of the weak problem. First, there is a sign error in the boundary term. If we derive a weak problem via integration by parts, we get the following: \int_\Omega -\nabla\cdot(k\nabla u)v\,d\mathbf{x} = \int_\Omega fv\,d\mathbf{x} \Rightarrow \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} ~-~\int_{\partial\Omega}(k\nabla u\cdot\mathbf{n})v\,d\mathbf{s}~=~\int_\Omega f v\,d\mathbf{x} \Rightarrow \int_\Omega k\nabla u\cdot\nabla v\,d\mathb…