# Problem with crank nicolson for coupled non linear pde

hi at all,

I started with FEA last December and my problem (thermoelectrcity) is not the best way to start with finite elements. But with the fenics and dolfin documentation and with some examples from the web I was able to solve the steady state problem.

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.

I think that the error is in the weak form and/or in the while loop.

But due to the lack of knowledge I find no way to solve my problem with the help of other documented examples.

I hope, somebody has time to look over my code and give me a right hint.

I had no idea how to reduce the problem into few lines.

I'm sorry. But maybe the problem is interesting enough to be used as a demo file. :)

best regards

Michael

Below the code:

import numpy

from dolfin import *

# Class representing the intial conditions

class InitialConditio

def __init_

def eval(self,vec,x):

def value_shape(self):

return (2,)

# Class for interfacing with the Newton solver

class ThermoelEquatio

def __init__(self, a, L,bc):

self.L = L

self.a = a

self.bc=bcs

def F(self, b, x):

for condition in self.bc:

def J(self, A, x):

for condition in self.bc:

# Model parameters

dt =5.0e-02 # time step

theta = 1# time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

Tinni = 300. # initial temperature

Vinni = 0.0 # initial electr. potential

# Form compiler options

parameters[

parameters[

parameters[

# Create mesh and define function spaces

# ===============

# Mesh of simple 2 dim. pellet

# ( 0.1 cm X 0.6 cm)

# ===============

mesh = RectangleMesh(

#######

#Subdomains

#three domains with different temperature depend. material properties

#######

class OmegaTop(

def inside(self, x, on_boundary):

return (between(x[1], (0.59, 0.6)) and between(x[0], (0.0, 0.1)))

class OmegaMiddle(

def inside(self, x, on_boundary):

return (between(x[1], (0.01, 0.59)) and between(x[0], (0.0, 0.1)))

class OmegaBottom(

def inside(self, x, on_boundary):

return (between(x[1], (0.0, 0.01)) and between(x[0], (0.0, 0.1)))

# Initialize sub-domain instances

omegatop = OmegaTop()

omegamiddle = OmegaMiddle()

omegabottom = OmegaBottom()

# Initialize mesh function for interior domains

subdomains = CellFunction(

subdomains.

omegatop.

omegamiddle.

omegabottom.

#plot(subdomain

#######

# boundary conditions

#######

class GammaBottom(

def inside(self, x, on_boundary):

return near(x[1], 0.0)

class GammaTop(

def inside(self, x, on_boundary):

return near(x[1], 0.6)

gammabottom =GammaBottom()

gammatop = GammaTop()

# Initialize mesh function for boundary domains

boundaries = FacetFunction(

boundaries.

gammabottom.

gammatop.

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)

te, vo = 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 = InitialConditio

u.interpolate(

u0.interpolate(

#I do not know if I'm right

temp_mid = (1.0-theta)

volt_mid = (1.0-theta)

#######

#Temperatur dependend material properties

#Coefficients of the polynomials A * temp² + B* temp + C

#######

#thermal conductivity

hcond1= Function(V1)

hcond2= Function(V1)

hcond3 = Function(V1)

#Seebeck coeff.

seeb1= Function(V1)

seeb2= Function(V1)

seeb3 = Function(V1)

#electr. conductivity

elcond1= Function(V1)

elcond2= Function(V1)

elcond3 = Function(V1)

k1_values = [0., 1.68181805028e-

k2_values = [0., -1.11448469117e

k3_values = [4., 0.0189634948598

a1_values = [0., -5.5012121918e-

a2_values = [0., 6.92678311478e-07 ,0.]

a3_values = [6.5e-6, -4.45128316556e

c1_values = [ 0.,0.0009090909

c2_values = [ 0.,-1.71636366069, 0.]

c3_values = [ 5.9e6,1603.

help = numpy.asarray(

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

hcond2.vector()[:] = numpy.choose(help, k2_values)

hcond3.vector()[:] = numpy.choose(help, k3_values)

elcond1.vector()[:] = numpy.choose(help, c1_values)

elcond2.vector()[:] = numpy.choose(help, c2_values)

elcond3.vector()[:] = numpy.choose(help, c3_values)

seeb1.vector()[:] = numpy.choose(help, a1_values)

seeb2.vector()[:] = numpy.choose(help, a2_values)

seeb3.vector()[:] = numpy.choose(help, a3_values)

#######

rho=7.8 #Density g/cm³

cp=0.28 #heat capacity J/gK

Tc = Constant(300.15) # 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. +800.*sin(

#V1 = Constant(0.2) # el. potential (top boundary)

#### Define boundary conditions

bc0 = DirichletBC(

bc1 = DirichletBC(

bc2 = DirichletBC(

#addionally

#bc4 = DirichletBC(

#and/ or Neumann:

n = FacetNormal(mesh)

qc = as_vector(

jc = as_vector(

g_T = dot(n,qc)

g_R = dot(n,jc)

bcs = [bc0,bc1,

ds = Measure(

#electric current flux per unit area A/cm²:

j = -(elcond1*

- (seeb1*

elcond2*

#thermal flux per unit area W/cm²:

q = -((hcond1*

+ (seeb1*

#electric field V/cm

E = -nabla_

# generated heat W/cm³ (inner source)

Q = dot(j,E)

# Weak statement of the equations

#VOLT: -inner(

#TEMP: -inner(

L0 = -inner(

L1 = -(rho*cp)

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.

solver.

# Output file

#fileU=

#fileT=

#######

# Step in time

t = 0.0

T = 20*dt

while (t < T):

print "======

print "Solving at t= "+str(t)+"s, tmax="+str(T)+\

"s, dt="+str(dt)+"s."

print "======

t += dt

u0.vector()[:] = u.vector()

solver.

#fileT << (u.split()[0], t)

#fileU << (u.split()[1], t)

plot(u.

plot(u.

interactive()

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- M Waste

- Solved:
- 2013-04-22

- Last query:
- 2013-04-22

- Last reply:
- 2013-04-19

Johan Hake (johan-hake) said : | #1 |

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.

Myles English (mylesenglish) said : | #2 |

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

M Waste (michael-waste) said : | #3 |

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

Jan Blechta (blechta) said : | #4 |

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

M Waste (michael-waste) said : | #5 |

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

Jan Blechta (blechta) said : | #6 |

On Fri, 08 Mar 2013 14:30:57 -0000

M Waste <email address hidden> wrote:

> Question #222911 on DOLFIN changed:

> https:/

>

> 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(

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

>

M Waste (michael-waste) said : | #7 |

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

Jan Blechta (blechta) said : | #8 |

On Fri, 08 Mar 2013 19:50:55 -0000

M Waste <email address hidden> wrote:

> Question #222911 on DOLFIN changed:

> https:/

>

> 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

>

M Waste (michael-waste) said : | #9 |

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

seebeck=

plot (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... :-)

Jan Blechta (blechta) said : | #10 |

On Sun, 10 Mar 2013 09:45:52 -0000

M Waste <email address hidden> wrote:

> Question #222911 on DOLFIN changed:

> https:/

>

> M Waste posted a new comment:

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

>

> seebeck=

> plot (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://

http://

?

>

> Therefore I have to read more special literature or how does one say,

> a cobbler should stick to his last... :-)

>

M Waste (michael-waste) said : | #11 |

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(

hcond1.

is this correct?

>Isn't this output what expected

>http://

>http://

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

ramzy daou (ramzy-daou) said : | #12 |

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(

L1 = -(rho*cp)

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

best regards,

Ramzy

M Waste (michael-waste) said : | #13 |

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

M Waste (michael-waste) said : | #14 |

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

ramzy daou (ramzy-daou) said : | #15 |

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 InitialConditio

def __init_

def eval(self,vec,x):

def value_shape(self):

return (2,)

# Class for interfacing with the Newton solver

class ThermoelEquatio

def __init__(self, a, L,bc):

self.L = L

self.a = a

self.bc=bcs

def F(self, b, x):

for condition in self.bc:

def J(self, A, x):

for condition in self.bc:

# Model parameters

dt =5.0e-03 # time step

Tinni = 300. # initial temperature

Vinni = 0.0 # initial electr. potential

# Form compiler options

parameters[

parameters[

parameters[

# Create mesh and define function spaces

# ===============

# Mesh of simple 2 dim. pellet

# ( 0.1 cm X 0.6 cm)

# ===============

mesh = RectangleMesh(

#######

# boundary conditions

#######

class GammaBottom(

def inside(self, x, on_boundary):

return near(x[1], 0.0)

class GammaTop(

def inside(self, x, on_boundary):

return near(x[1], 0.6)

gammabottom =GammaBottom()

gammatop = GammaTop()

# Initialize mesh function for boundary domains

boundaries = FacetFunction(

boundaries.

gammabottom.

gammatop.

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 = InitialConditio

u.interpolate(

u0.interpolate(

#######

#Temperatur dependend material properties

#Coefficients of the polynomials A * temp² + B* temp + C

#######

seebeck0 = Constant(

seebeck1 = Constant(

seebeck2 = Constant(

sigma0 = Constant(1600.0)

sigma1 = Constant(

sigma2 = Constant(

kappa0 = Constant(

kappa1 = Constant(

kappa2 = Constant(

#######

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(

#### Define boundary conditions

bc0 = DirichletBC(

bc1 = DirichletBC(

bc2 = DirichletBC(

bcs = [bc0,bc1,

ds = Measure(

# 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*

+ dt*s*seebeck*

- dt*inner(

- dt*inner(

- s*cp*rho*temp*dx + s*cp*rho*

L1 = - inner(nabla_

- inner(nabla_

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.

solver.

# Output file

#fileU=

#fileT=

#######

# 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.

#fileT << (u.split()[0], t)

#fileU << (u.split()[1], t)

plot(u.

plot(u.

interactive()

M Waste (michael-waste) said : | #16 |

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

M Waste (michael-waste) said : | #17 |

sorry I meant:

...

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

Best regards

Michael

ramzy daou (ramzy-daou) said : | #18 |

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

Jan Blechta (blechta) said : | #19 |

On Fri, 19 Apr 2013 09:07:08 -0000

ramzy daou <email address hidden> wrote:

> Question #222911 on DOLFIN changed:

> https:/

>

> 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

>

M Waste (michael-waste) said : | #20 |

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

Jan Blechta (blechta) said : | #21 |

On Mon, 22 Apr 2013 09:35:57 -0000

M Waste <email address hidden> wrote:

> Question #222911 on DOLFIN changed:

> https:/

>

> 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

>

ramzy daou (ramzy-daou) said : | #22 |

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

M Waste (michael-waste) said : | #24 |

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

ramzy daou (ramzy-daou) said : | #25 |

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://

This has already been modelled, for example in Comsol by Jaegle:

www.comsol.

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

Jan Blechta (blechta) said : | #26 |

On Fri, 03 May 2013 08:31:19 -0000

ramzy daou <email address hidden> wrote:

> Question #222911 on DOLFIN changed:

> https:/

>

> 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://

> This has already been modelled, for example in Comsol by Jaegle:

> www.comsol.

>

> 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://

http://

Jan