[Getdp] ... Formulation question ...

mkoch at gvtc.com mkoch at gvtc.com
Thu Mar 8 16:09:15 CET 2007


Hi Olivier, Hi Bernhard,

I was able to make the solution "move" by introducing a second  
Quantity for the test function. So note the use of Lvl AND LvlTst below.

Formulation {
   {Name FrmLvl; Type FemEquation;
     Quantity {
       {Name Lvl; Type Local; NameOfSpace SpcLvl;}
       {Name LvlTst; Type Local; NameOfSpace SpcLvl;}
     }
     Equation {
       Galerkin {Dt[Dof{Lvl},{LvlTst}];In Vlm;Integration Int;Jacobian Jac;}
       Galerkin {[Spd[]*Norm[{d Lvl}],{LvlTst}];In Vlm;Integration  
Int;Jacobian Jac;}
//      Galerkin {[Vel[]*Dof{d Lvl},{LvlTst}];In Vlm;Integration  
Int;Jacobian Jac;}
     }
   }
}

Unfortunately, the "solution" to this now is completely unphysical  
with huge numbers, like 10^30 or thereabouts! This is a real puzzler!  
Even more so, because when you comment out the "Spd" Galerkin term and  
uncomment the "Vel" Galerkin term, it works nicely, except for the  
aforementioned jaggedness.

Regards,

Matt Koch


----- Message from mkoch at gvtc.com ---------
     Date: Thu,  8 Mar 2007 08:17:45 -0600
     From: mkoch at gvtc.com
Reply-To: mkoch at gvtc.com
  Subject: Re: [Getdp] ... Formulation question ...
       To: getdp at geuz.org


> Hi Bernhard,
>
> thanks for your support. I think I do exactly as you said in my
> Constraint{}. I use the following line:
>
> {Region VlmPlt;Type Init;Value X[];}
>
> If I am not too mistaken, in 2D, this should generate a plane that
> inclines linearly along the x-axis, and the "solution" actually bears
> that out. Also:
>
> d/dt(Phi) - v * Grad(Phi) = 0
>
> works perfectly fine, even though the solution becomes very jagged at
> the edges after a while, but:
>
> d/dt(Phi) - s * Norm(Grad(Phi)) = 0
>
> does not seem to advance in time. And the Galerkin term containing the
> velocity vector v vs. the one containing the speed scalar s is
> literally the only difference between the two! Well that, and using
> Dof{d Phi} in the velocity term and {d Phi} in the speed term.
>
> Therefore, I believe this has something to do with the usage or
> non-usage of "Dof". I think Olivier Castany is checking out my files,
> so we'll see what he might find.
>
> Best Regards,
>
> Matt Koch
>
>
>
> ----- Message from Bernhard.Kubicek at arsenal.ac.at ---------
>      Date: Thu, 8 Mar 2007 10:04:35 +0100
>      From: Kubicek Bernhard <Bernhard.Kubicek at arsenal.ac.at>
> Reply-To: Kubicek Bernhard <Bernhard.Kubicek at arsenal.ac.at>
>   Subject: Re: [Getdp] ... Formulation question ...
>        To: mkoch at gvtc.com, getdp at geuz.org
>
>
>> Hello!
>> Maybe the reason is in your formula:
>> d/dt(Phi)= Spd| grad phi|
>> if phi is initialized with zero, it will always be zero (apart from
>> assigned boundary conditions).
>> Maybe you can initalize phi with something usefull,e.g. using a
>> coordinate dependent function in the constraint as value. Watch out
>> for the difference betweeen "assign" and "init" constraints.
>> Another possibility, would be to calculate an initial solution for
>> phi by its own formulation, and then use a chained resolution to
>> transfer this initial field to your already existing transient
>> solution.
>> nice greetings,
>>  bernhard
>>
>>
>> -----Ursprüngliche Nachricht-----
>> Von: getdp-bounces at geuz.org [mailto:getdp-bounces at geuz.org] Im
>> Auftrag von mkoch at gvtc.com
>> Gesendet: Mittwoch, 7. März 2007 18:33
>> An: getdp at geuz.org
>> Betreff: Re: [Getdp] ... Formulation question ...
>>
>>
>> Hello Olivier,
>>
>> so I tried exactly what you suggested below (I think). And I tried it
>> with the TimeLoopTheta scheme set to 0.0, 0.5 and 1.0. In all cases,
>> GetDP no longer complains about any errors, which is nice. BUT, in
>> neither case does it appear as if the second Galerkin term is anything
>> but zero. I am certain that Spd[] is not zero, so I suspect Norm[{d
>> Phi}] somehow remains zero?
>>
>> I even tried Sqrt[{d Phi}*{d Phi}] instead of Norm[{d Phi}], but that
>> does not work either. How is that possible? In the general heat
>> conduction equation posted on the Wiki, there is a radiation term, and
>> it uses something equivalent to {Phi}^4, and that seems to work
>> (although, it also uses source and boundary terms, which are absent
>> from my equations).
>>
>> How come {Phi}^4 turns into something non-zero after a while, but
>> Norm[{d Phi}] seemingly does not? I mean, I even know that Phi is
>> initialized to something that does not create a zero Norm[{d Phi}].
>> So, from the start, Norm[{d Phi}] should be noon-zero and thus drive
>> Dt[Dof{Phi},{Phi}] to something other than zero for the other time
>> steps as well!? Instead, I get a solution that is equivalent to
>> d/dt(Phi) = 0, rather than the desired d/dt(Phi) + Spd*Norm(Grad(Phi))
>> = 0.
>>
>> Any thoughts?
>>
>> Regards,
>>
>> Matt
>>
>> ----- Message from castany at quatramaran.ens.fr ---------
>>      Date: Sun, 4 Mar 2007 21:37:55 +0100
>>      From: Olivier Castany <castany at quatramaran.ens.fr>
>> Reply-To: Olivier Castany <castany at quatramaran.ens.fr>
>>   Subject: Re: [Getdp] ... Formulation question ...
>>        To: getdp at geuz.org
>>
>>
>>> Hello,
>>>
>>> now that you give a more precise desciption of your problem, I
>>> understand that you want to solve a dynamic equation.
>>>
>>> I can only be of little help since I do not understand the part b) of
>>> the algorithm you want to program.
>>>
>>> I can just say that it is possible to program the equation :
>>>
>>>> a) d/dt(Phi) + Spd*Norm(Grad(Phi)) = 0
>>>
>>> the formulation would be :
>>>
>>> Galerkin { Dt [ Dof{Phi} , {Phi} ] ; ... }
>>> Galerkin { [ Spd[] * Norm[{d Phi}]  , {Phi} ] ; ... }
>>>
>>> and the resolution with "TimeLoopTheta" will be an explicit Euler
>>> algorithm.
>>>
>>> As I don't really understand your algorithm with b), I can't say if it
>>> is possible to do it with GetDP. My feeling is that it is not
>>> currently possible without modifying the program.
>>>
>>> Regards,
>>>
>>> O.C.
>>>
>>> PS : if you manage to do something nice with the level set equations,
>>> maybe you can post it later or put it on the wiki ?
>>> _______________________________________________
>>> getdp mailing list
>>> getdp at geuz.org
>>> http://www.geuz.org/mailman/listinfo/getdp
>>>
>>
>>
>> ----- End message from castany at quatramaran.ens.fr -----
>>
>>
>>
>>
>>
>> _______________________________________________
>> getdp mailing list
>> getdp at geuz.org
>> http://www.geuz.org/mailman/listinfo/getdp
>>
>> _______________________________________________
>> getdp mailing list
>> getdp at geuz.org
>> http://www.geuz.org/mailman/listinfo/getdp
>>
>
>
> ----- End message from Bernhard.Kubicek at arsenal.ac.at -----
>
>
>
>
>
> _______________________________________________
> getdp mailing list
> getdp at geuz.org
> http://www.geuz.org/mailman/listinfo/getdp
>


----- End message from mkoch at gvtc.com -----