Subject: [Unicorn-dev] What's everyone working on?

From: Peter Brune <[email-masked]>
Date: Tue, 11 Aug 2009 11:45:14 -0500

I guess I'm being too polite here and not pushing things I see get edited; I
think we should break up the lib dir and just fix it; I'm going to start
working on the elasticity stuff -- it has forms to update and I can do this
quickly.

Anything else?

- Peter

---

From: Niclas Jansson <[email-masked]>
Date: Tue, 11 Aug 2009 20:16:12 +0200

> working on the elasticity stuff -- it has forms to update and I can do this

Isn't only TimeDependent left to modify before everything compiles? 
Once we have the icns solver working, I'll start debugging slip.

Niclas

---

From: Peter Brune <[email-masked]>
Date: Tue, 11 Aug 2009 13:23:56 -0500

> Once we have the icns solver working, I'll start debugging slip.

There's all sorts of nastiness with regard to mesh quality; I've started
just moving things out of the lib directory in my local build.  The new
organization of icns is something we might want to all discuss tomorrow
too.  I've been looking at it and have slowly started to piece it together.
I'm assuming that NewSlipBC, NewDirichletBC, and Project are redundant at
this point?

---

From: Johan Jansson <[email-masked]>
Date: Wed, 12 Aug 2009 02:27:12 +0200

Peter Brune wrote:
>
>
> On Tue, Aug 11, 2009 at 1:16 PM, Niclas Jansson <[email-masked] 
> <mailto:[email-masked]>> wrote:
>
>     Peter Brune <[email-masked] <mailto:[email-masked]>> writes:
>
> There's all sorts of nastiness with regard to mesh quality; I've 
> started just moving things out of the lib directory in my local 
> build.  The new organization of icns is something we might want to all 
> discuss tomorrow too.  I've been looking at it and have slowly started 
> to piece it together.  I'm assuming that NewSlipBC, NewDirichletBC, 
> and Project are redundant at this point?
Hi,

What's the nastiness with regard to mesh quality? That it doesn't 
compile, or something else? NewSlipBC is needed for a Newton-type method 
where the Jacobian may not always be updated. This can wait until 
everything else is working.

I saw you defined stabilization functions as C++ expressions. Isn't it 
better to aim to put all of that into UFL, something like:

v      = TestFunction(vector)      # test basis function
U      = Function(vector)     # trial basis function

k  = Constant(tetrahedron) # time step
h = Constant(tetrahedron) # cell size

delta1 = 4.0 * (0.5 / sqrt(1.0 / (k**2)) + sqrt(inner(U, U))/h**2)

or?

  Johan

---

From: Peter Brune <[email-masked]>
Date: Wed, 12 Aug 2009 01:28:58 -0500

> What's the nastiness with regard to mesh quality? That it doesn't

It doesn't compile and has terrible, terrible programming practices that
make it seem not worth making compile. :)

> delta1 = 4.0 * (0.5 / sqrt(1.0 / (k**2)) + sqrt(inner(U, U))/h**2)

Isn't this wrong?  This is a global quantity, whereas the stabilization you
see in the ICNS solver changes locally with h, nu, and depends on the global
norm of U.  So, it's a combination of a global and local evaluation and I
saw it best done as a function.

- Peter

---

From: Johan Jansson <[email-masked]>
Date: Wed, 12 Aug 2009 10:12:07 +0200

Peter Brune wrote:
> I guess I'm being too polite here and not pushing things I see get 
> edited; I think we should break up the lib dir and just fix it; I'm 
> going to start working on the elasticity stuff -- it has forms to 
> update and I can do this quickly.
>
> Anything else?
>
> - Peter
> ------------------------------------------------------------------------
>
> _______________________________________________
> Unicorn-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/unicorn-dev
>   
Hi,

Please don't hesitate to take a look at re-organizing the directory 
structure if you want. Don't worry too much about being polite at this 
stage :). most things are up in the air now so now is the time to do big 
changes.

  Johan

---

From: Johan Jansson <[email-masked]>
Date: Wed, 12 Aug 2009 10:38:29 +0200

Peter Brune wrote:
> On Tue, Aug 11, 2009 at 7:27 PM, Johan Jansson <[email-masked] 
> <mailto:[email-masked]>> wrote:
>
>     Peter Brune wrote:
>     >
>     >
>     > On Tue, Aug 11, 2009 at 1:16 PM, Niclas Jansson
>     <[email-masked] <mailto:[email-masked]>
>     > <mailto:[email-masked] <mailto:[email-masked]>>> wrote:
>     >
>     >     Peter Brune <[email-masked] <mailto:[email-masked]>
>     <mailto:[email-masked] <mailto:[email-masked]>>> writes:
>     >
>     >     > think we should break up the lib dir and just fix it; I'm
>     going
>     >     to start
>     >     > working on the elasticity stuff -- it has forms to update
>     and I
>     >     can do this
>     >     > quickly.
>     >     >
>     >     > Anything else?
>     >     >
>     >
>     >
>     >
>     >     > - Peter
>     >
> It doesn't compile and has terrible, terrible programming practices 
> that make it seem not worth making compile. :) 
I don't think you are talking about the MeshQuality class, or? Perhaps 
you are talking about the ElasticSmoother class? This is quite messy, 
and probably a better idea to rewrite before porting to the new 
interface. The MeshQuality class is very clean though and already 
compiles as far as I can see.

> Isn't this wrong?  This is a global quantity, whereas the 
> stabilization you see in the ICNS solver changes locally with h, nu, 
> and depends on the global norm of U.  So, it's a combination of a 
> global and local evaluation and I saw it best done as a function.
>
> - Peter
>
It's supposed to be the (local) Euclidean norm of U, not the global norm 
(U is the local convective velocity).

You're right that h is supposed to be defined as a Function though, so:

k  = Constant(tetrahedron) # time step
h = Function(constant_scalar) # cell size

The thing we're missing here is the switch on h and nu:

if  h/nu > 1 or ny < 10^-10
  d1 = C1 * ( 0.5 / sqrt( 1/k^2 + |U|^2/h^2 ) )
else
  d1 = C1 * h^2

UFL has support for conditionals, so we could implement it directly. 
However, I don't think it actually works in FFC yet.

For the high Reynolds number applications which is our main focus, then 
this will always evaluate to the first case, so as a first 
implementation we can just skip the second case. Then we can think how 
to implement this conditional, perhaps this could also be done without 
the conditional with some approximation of a step function.

But let's try to do this in the form? This would reduce the code size 
and increase maintainability.

  Johan

---

From: Anders Logg <[email-masked]>
Date: Wed, 12 Aug 2009 11:51:32 +0200

> However, I don't think it actually works in FFC yet.

No, it doesn't but it would be good to have to support things like
this. You could get around it by adding some weight function which is
either 0 or 1 depending on the value of h and nu, and then weight the
two terms with w and (1 - w).

--
Anders