Problem with crank nicolson for coupled non linear pde

Not many on this forum will debug your full code. Cut down the problem
until it is small and comprehensible and then it might be that you find
the answer your self. If not post it here again.

Johan

On 02/26/2013 10:15 PM, M Waste wrote:
> Now I have no idea how I should implement the CN-method to solve the coupled problem in transient mode.
> I used the cahn hilliard example for my problem but I do not get a correct time progress.
> The solution is only a switching between the initial condition and the final time step.
> Regardless of the time step size.

Hi Michael,

I don't know what difference it will make but I think the line
self.bc=bcs should be self.bc=bc

Myles

hello Johan, hello Myles,

thank you for answering.
Of course a small sample code would be helpful.
And you are right no one has time to dive into my problem in detail.
Unfortunately I'm not able to cut down the problem. I tried it before I have posted the problem .
But I got always a non running code.
I'm a noob and I was glad to see a reasonable result in steady state case despite my cumbersome programming.

This code run without an error message. But not correct.
That's why I though that an experienced user will see the mistake with low effort.

@ Myles: Thank you for the hint. I changed it, but it has no effect.

I found a non linear coupled example with an "update function" which is called inside the while loop. I'll try to adapt it for my problem. Maybe that works.

But I'm open for further help.

Michael

hcond1.vector()[:] = numpy.choose(help, k1_values)

You assume dofs ordering. Dofs may not be no more ordered as some mesh entity in DOLFIN 1.1.0+. Try looking to some recent answers regarding 'dofs ordering'. There are plenty of them.

Jan

Hi Jan,

thank you for your reply.
I'm not sure what you mean.
In steady state case it was not a problem with this code. The solution was the same as a literature example.

However I have a problem with the time depended weak form of the thermoelectric pde.
I'm quite sure that I missed something.

I need more practice to understand the mathematics behind and also the FEA methods.

I'm so sorry for bothering you with my problem.

Michael

On Fri, 08 Mar 2013 14:30:57 -0000
M Waste <email address hidden> wrote:
> Question #222911 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/222911
>
> Status: Answered => Open
>
> M Waste is still having a problem:
> Hi Jan,
>
> thank you for your reply.
> I'm not sure what you mean.

hcond1.vector()[:] is indexed by indices called DOF (degree of freedom)
numbers but

help = numpy.asarray(subdomains.array(), dtype=numpy.int32) is indexed
by cell indices

There is no correspondence between these indices since DOLFIN 1.1.0.
Check questions 222203 and 221862.

> In steady state case it was not a problem with this code. The
> solution was the same as a literature example.

Maybe you are using DOLFIN 1.0.0. Then my comment does not apply.

>
> However I have a problem with the time depended weak form of the
> thermoelectric pde. I'm quite sure that I missed something.
>
> I need more practice to understand the mathematics behind and also the
> FEA methods.
>
> I'm so sorry for bothering you with my problem.
>
> Michael
>

hello Jan,

I checked the installed version via python -c"import dolfin;print dolfin.__version__"
and I got 1.1.0.

I'm really confused. I used this code line, because it was for me the only "simple" way to load different subdomains with temperature depended material properties.
It is a stupid questions, but why ist this not a problem in steady state? Or maybe it is and I don't see it.

Michael

On Fri, 08 Mar 2013 19:50:55 -0000
M Waste <email address hidden> wrote:
> Question #222911 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/222911
>
> M Waste posted a new comment:
> hello Jan,
>
> I checked the installed version via python -c"import dolfin;print
> dolfin.__version__" and I got 1.1.0.
>
> I'm really confused. I used this code line, because it was for me the
> only "simple" way to load different subdomains with temperature
> depended material properties. It is a stupid questions, but why ist
> this not a problem in steady state? Or maybe it is and I don't see
> it.

Try plotting these material coefficients if you can see what's
expected. Speaking generally - you should test every small part
of code if doing what expected.

Jan

>
>
> Michael
>

I checked the material coeff. with the following code line.

seebeck=project(seeb1*temp**2+seeb2*temp+seeb3,FunctionSpace(mesh,"CG",1))
plot (seebeck,scale=0.,title='seebeck')
and also for heat conductivity and el. conductivity. The plots look as expected. Also the "boundary" between the subdomains show no deviations or discontinuities.

However I'm not able to check the time dependend part, because I'm not (yet) familiar with weak formulations.
I think the Crank Nicolson part in my code is more or less nonsense.

Therefore I have to read more special literature or how does one say, a cobbler should stick to his last... :-)

On Sun, 10 Mar 2013 09:45:52 -0000
M Waste <email address hidden> wrote:
> Question #222911 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/222911
>
> M Waste posted a new comment:
> I checked the material coeff. with the following code line.
>
> seebeck=project(seeb1*temp**2+seeb2*temp+seeb3,FunctionSpace(mesh,"CG",1))
> plot (seebeck,scale=0.,title='seebeck')
> and also for heat conductivity and el. conductivity. The plots look
> as expected. Also the "boundary" between the subdomains show no
> deviations or discontinuities.

Ok, I see. But it is working "by accident". One should no more make
these assumptions. For exmaple this coincidence is no more valid for
CG1 dofs and vertices. But probably it is still working for DG0 dofs
and cells.

>
> However I'm not able to check the time dependend part, because I'm
> not (yet) familiar with weak formulations. I think the Crank
> Nicolson part in my code is more or less nonsense.

Isn't this output what expected
http://artax.karlin.mff.cuni.cz/~blecj6am/222911/dolfin_plot_0.png
http://artax.karlin.mff.cuni.cz/~blecj6am/222911/dolfin_plot_1.png
?

>
> Therefore I have to read more special literature or how does one say,
> a cobbler should stick to his last... :-)
>

Thank you very much for your support!

>Ok, I see. But it is working "by accident". One should no more make
>these assumptions. For exmaple this coincidence is no more valid for
>CG1 dofs and vertices. But probably it is still working for DG0 dofs
>and cells.

I see. I'll try to change this im my code.
as seen in question 221862
for cell in cells(mesh):
   subdomain_no = subdomains[cell]
   dof_ind = V1.dofmap().cell_dofs(cell.index())[0]
   hcond1.vector()[dof_ind] = k1_values[subdomain_no]

is this correct?

>Isn't this output what expected
>http://artax.karlin.mff.cuni.cz/...
>http://artax.karlin.mff.cuni.cz/....

Thank you for this print out. For the temperature I got the same plot. but for the el. potential is the value range between 0 and -0.100.

How did you get this solution?

Michael

Hello Michael,

I am very interested in implementing the thermoelectric problem in Fenics as well. However, it seems quite tricky, given the non-linearity. Did you make any further progress?

As regards the time-dependent part of the problem, I think you are missing a couple of terms. I believe temp_mid and temp_prev should enter into L1 like so:

L0 = -inner(nabla_grad(vo),j)*dx + g_R*vo*ds(2)
L1 = -(rho*cp)*te*temp_mid*dx +(rho*cp)*te*temp_prev*dx - dt*inner(nabla_grad(te),q)*dx +dt*Q*te*dx #+ te*g_T*ds(1)

But I confess I am also a beginner with FE and Fenics too.

best regards,
Ramzy

Hello Ramzy,

unfortunately I hadn't enough time in the last weeks.
Thank you for the hint.
I going to implement the missing terms and let you know how it works.

best regards,

Michael

Hello Ramzy,

short summary of my trial.
with this L1 term I haven't reached convergence.
I switched the signs and then I reached convergence. But the code shows the same behaviour.
There is no progress in time.
Unfortunately I do not see the mistake.
If you are able to solve the issue please let me know about your progress.

Best regards,

Michael

Hi Michael,

I had another look today. I think you also had a problem with your initial conditions, Uini and Tini are in the wrong place in the class definition.

I have tried expanding completely the equations so that they depend only on T and V. Below some code where I have considered only a single domain for simplicity. I have tried changing the time step and the results seem sensible to me, but let me know if it matches your textbook.

For Dirichlet boundary conditions only (or zero-flux Neumann), this seems to work OK. I am now not sure how to properly include Neumann conditions, e.g. for an electrical current flowing through the sample, given that current is not well defined any more. Any ideas?

best regards,
Ramzy

# -*- coding: utf-8 -*-
"""
Michael Waste's code for thermoelectric FE in fenics
Modified by Ramzy Daou 15 April 2013
"""

import numpy
from dolfin import *

# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self,Tini,Uini):
        self.Tini=Tini
        self.Uini=Uini
    def eval(self,vec,x):
        vec[0]=self.Tini#T ( 0 ) = T i n i
        vec[1]=self.Uini#u( 0 ) = 0.
    def value_shape(self):
        return (2,)
# Class for interfacing with the Newton solver
class ThermoelEquation(NonlinearProblem):
    def __init__(self, a, L,bc):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bc=bcs
        self.reset_sparsity = True
    def F(self, b, x):
        assemble(self.L, tensor=b)
        for condition in self.bc:
            condition.apply(b,x)
    def J(self, A, x):
        assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity)
        for condition in self.bc:
            condition.apply(A)
        self.reset_sparsity = False

# Model parameters

dt =5.0e-03 # time step
Tinni = 300. # initial temperature
Vinni = 0.0 # initial electr. potential
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"

# Create mesh and define function spaces
# ===============
# Mesh of simple 2 dim. pellet
# ( 0.1 cm X 0.6 cm)
# ===============
mesh = RectangleMesh(0,0,0.1,0.6,10,60)

#####################################
# boundary conditions
#####################################
class GammaBottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)
class GammaTop(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.6)

gammabottom =GammaBottom()
gammatop = GammaTop()
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh,2)
boundaries.set_all(0)
gammabottom.mark(boundaries, 1)
gammatop.mark(boundaries, 2)

V1 = FunctionSpace(mesh, 'DG', 0)# temperatur depend. properties
P1 = FunctionSpace(mesh, "Lagrange", 1)#temperature
P2 = FunctionSpace(mesh, "Lagrange", 1)#electric potential

ME = P1*P2

# Define trial and test functions
du = TrialFunction(ME)
s, r = TestFunctions(ME)

# Define functions
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step

# Split mixed functions
dc, dp = split(du)
temp, volt = split(u)
temp_prev, volt_prev = split(u0)#previous step

# Create intial conditions and interpolate
u_init = InitialConditions(Tinni,Vinni)
u.interpolate(u_init)
u0.interpolate(u_init)

#####################################################################
#Temperatur dependend material properties
#Coefficients of the polynomials A * temp² + B* temp + C
#####################################################################

seebeck0 = Constant(-4.45128316556e-05)
seebeck1 = Constant(6.92678311478e-07 )
seebeck2 = Constant(-5.5012121918e-10)

sigma0 = Constant(1600.0)
sigma1 = Constant(-1.71636366069)
sigma2 = Constant(0.000909090929458)

kappa0 = Constant(0.0189634948598)
kappa1 = Constant(-1.11448469117e-05)
kappa2 = Constant(1.68181805028e-08)

######################################################################
rho=7.8 #Density g/cm³
cp=0.28 #heat capacity J/gK
Tc = Constant(300.00) # cold side temperature (bottom boundary)
V0 = Constant(0.0) # el. potential (bottom boundary)
Th = Constant(873.15) # hot side temperature (top boundary)
# Th could alsovary in time like:
#Th = Expression('400. +100.*sin(omega*t)', omega=omega, t=0.0)
#### Define boundary conditions
bc0 = DirichletBC(ME.sub(0), Th, boundaries, 2)
bc1 = DirichletBC(ME.sub(1), V0, boundaries, 1)
bc2 = DirichletBC(ME.sub(0), Tc, boundaries, 1)

bcs = [bc0,bc1,bc2]#,bc3,bc2,bc4bc0

ds = Measure("ds")[boundaries]

# thermoelectric power
seebeck = seebeck0 + seebeck1*temp +seebeck2*temp*temp
#thermal conductivity
kappa = kappa0 + kappa1*temp + kappa2*temp*temp
# electrical conductivity
sigma = sigma0 + sigma1*temp + sigma2*temp*temp

L0 = dt*s*sigma*inner(nabla_grad(volt),nabla_grad(volt))*dx\
     + dt*s*seebeck*sigma*inner(nabla_grad(volt),nabla_grad(temp))*dx\
     - dt*inner(nabla_grad(s),(kappa+seebeck*seebeck*sigma*temp)*nabla_grad(temp))*dx\
     - dt*inner(nabla_grad(s),(seebeck*sigma*temp)*nabla_grad(volt))*dx\
     - s*cp*rho*temp*dx + s*cp*rho*temp_prev*dx

L1 = - inner(nabla_grad(r),(seebeck*sigma)*nabla_grad(temp))*dx\
     - inner(nabla_grad(r),(sigma)*nabla_grad(volt))*dx

L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = ThermoelEquation(a, L,bcs)
solver = NewtonSolver("lu")
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
#fileU=File("speich/send_volt2dim.pvd")
#fileT=File("speich/send_temp2dim.pvd")
################################
# Step in time
t = 0.0
T = 100*dt
while (t < T):
    print "================================================================="
    print "Solving at t= "+str(t)+"s, tmax="+str(T)+\
    "s, dt="+str(dt)+"s."
    print "================================================================="
    t += dt
    V1.t = t
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    #fileT << (u.split()[0], t)
    #fileU << (u.split()[1], t)

plot(u.split()[0],scale=0.,title='temperature')
plot(u.split()[1],scale=0.,title='potential')

interactive()

hello Ramzy,

thanks a lot.
 I think your proposal is the right way.

Neumann: Did you try the following:
e.g.
jc = Constant(-20000.) #current density per square cm

and add this to the L1 term: ...+ dt*r*g_R*ds(2)

Best regards

Michael

sorry I meant:
...
and add this to the L1 term: ...+ dt*r*jc*ds(2)

Best regards

Michael

Hi Michael,

Thanks! It really is so simple. I was looking for a complicated solution. Well, I think we have a good working implementation of thermoelectric transport in Fenics. Can you think of any other features specific to the problem that we should include?

best regards,
Ramzy

On Fri, 19 Apr 2013 09:07:08 -0000
ramzy daou <email address hidden> wrote:
> Question #222911 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/222911
>
> ramzy daou posted a new comment:
> Hi Michael,
>
> Thanks! It really is so simple. I was looking for a complicated
> solution. Well, I think we have a good working implementation of
> thermoelectric transport in Fenics. Can you think of any other
> features specific to the problem that we should include?

Well, you could consider contributing a demo.

Jan

>
> best regards,
> Ramzy
>

First of all. Thank you both for your help.

At the moment I do not see any further specific features.
A coupling with the heat exchanger heat flow in order to simulate the whole integration is maybe possible but
for me to complex.

Demo: I'm mot the right person for that. My programming is to cumbersome as you have seen in the last month. :-)
My code was more or less useless without Jans hints and Ramzys time dependend code lines.
But generally I have no problem to contribute the above code as demo.

best regards,

Michael

On Mon, 22 Apr 2013 09:35:57 -0000
M Waste <email address hidden> wrote:
> Question #222911 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/222911
>
> Status: Answered => Solved
>
> M Waste confirmed that the question is solved:
> First of all. Thank you both for your help.
>
> At the moment I do not see any further specific features.
> A coupling with the heat exchanger heat flow in order to simulate the
> whole integration is maybe possible but for me to complex.
>
> Demo: I'm mot the right person for that. My programming is to
> cumbersome as you have seen in the last month. :-) My code was more
> or less useless without Jans hints and Ramzys time dependend code
> lines. But generally I have no problem to contribute the above code
> as demo.

Well, If you make it much simpler and shorter, to the most simpler
non-trivial problem of thermoelectric transport, it could be pushed
directly or I (or someone else...) could comb its hair if needed.

Jan

>
> best regards,
>
>
> Michael
>

Hello Jan,

I too am happy for my contribution to be included as an example. I think it is in the simplest form already that reflects a real-world example (single domain with Dirichlet boundary conditions, 1st order backward difference time dependence). I have tidied the code a little bit and also written a quick pdf with the source equations and variational form. I will email them to you.

best regards,
Ramzy

Hello Ramzy,

thank you for your summarization.
I tried it during the past week but good to see that I can stop LaTeXing :-)

Best regards

Michael

Hello Michael, Jan,

We could of course (as Jan has suggested) explore some examples that produce more interesting output. Here are a few suggestions:

Static problems:
1) The classic configuration for measuring the Seebeck effect in the steady state (with no electrical load)
Fixed j=0
applied q != 0
V = T = 0 at the sink (boundary 1)

2) The classic Peltier element: find the temperature rise (with no thermal load) for a given applied current
Fixed q = 0
applied j != 0
V = T = 0 at the sink

3) The loaded generator / heater / cooler
Include an electrical or thermal load in series with the thermoelectric element. How much power can be generated? How much heating/cooling? What is the device efficiency? A 'real-world' situation.

Dynamic problems:
4) Explore oscillatory currents or boundary conditions to look at the frequency response.

5) Transient pulse supercooling of a Peltier element as suggested here:
http://jap.aip.org/resource/1/japiau/v92/i3/p1564_s1
This has already been modelled, for example in Comsol by Jaegle:
www.comsol.com/papers/5256/download/Jaegle.pdf

6) Include charging effects in dielectric media by including the displacement field in the current equation.
This will make the time dependence more complicated.

I suppose 5) could produce some interesting output, while 3) has a lot of immediate applications.

best regards,
Ramzy

On Fri, 03 May 2013 08:31:19 -0000
ramzy daou <email address hidden> wrote:
> Question #222911 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/222911
>
> ramzy daou posted a new comment:
> Hello Michael, Jan,
>
> We could of course (as Jan has suggested) explore some examples that
> produce more interesting output. Here are a few suggestions:
>
> Static problems:
> 1) The classic configuration for measuring the Seebeck effect in the
> steady state (with no electrical load) Fixed j=0
> applied q != 0
> V = T = 0 at the sink (boundary 1)
>
> 2) The classic Peltier element: find the temperature rise (with no
> thermal load) for a given applied current Fixed q = 0
> applied j != 0
> V = T = 0 at the sink
>
> 3) The loaded generator / heater / cooler
> Include an electrical or thermal load in series with the
> thermoelectric element. How much power can be generated? How much
> heating/cooling? What is the device efficiency? A 'real-world'
> situation.
>
> Dynamic problems:
> 4) Explore oscillatory currents or boundary conditions to look at the
> frequency response.
>
> 5) Transient pulse supercooling of a Peltier element as suggested
> here: http://jap.aip.org/resource/1/japiau/v92/i3/p1564_s1
> This has already been modelled, for example in Comsol by Jaegle:
> www.comsol.com/papers/5256/download/Jaegle.pdf
>
> 6) Include charging effects in dielectric media by including the
> displacement field in the current equation. This will make the time
> dependence more complicated.
>
> I suppose 5) could produce some interesting output, while 3) has a lot
> of immediate applications.
>
> best regards,
> Ramzy
>

Well, this seems promising. My point is that code contains very much
boilerplate (many solution dependent coeeficients, complicated
form, ...) to produce non-interesting solution. NonlinearProblem
subclassing is already demonstrated in cahn-hilliard demo with much
more interesting output. Every non-linearity which is not needed to
produce non-trivial results could be canceled.

To make it more interesting you could also try to include some
non-trivial mesh from dolfin mesh repository, like dolphin:) Search
them here:
http://fenicsproject.org/download/data.html
http://fenicsproject.org/pub/data/meshes/

Jan