| Author | Message | ||
     Fred Rodriguez (plasmon) New member Username: plasmon Post Number: 1 Registered: 12-2007 |
Hi. Can anyone provide me with assistance? I am trying to understand why the following diffusion equation (see input file) works with the default FlexPDE BC (i.e., Natural(expr) = 0), with Natural BC provided by me but not with the Dirichlet BC. In the boundaries section, I've commented the lines that worked and left the Dirichlet BC uncommented. Whenever I run the code with the Dirichlet BC I get a flatline which is not to be expected. My end goal is to test the behavior of this Gaussian with different BCs because later I will expand this problem into a 3D problem and I am looking for the best way to define the BCs. Thanks! -------- Input File ------- TITLE 'Gaussian BC test' COORDINATES cartesian1 DEFINITIONS Lx=1e3 sigma =1 x0 = Lx/2 tmax=100 k=1 VARIABLES rho (threshold = 1e-2) INITIAL VALUES rho=exp(-(x)**2/(2*sigma^2))/(2*PI*sigma^2) EQUATIONS dt(rho)=k*dxx(rho) BOUNDARIES region 'domain' start (0) value(rho) = val(rho,0) line to (Lx) value(rho) = val(rho,lx) ! Dirichlet ! start (0) natural(rho) = -normal(0,k*dx(rho)) line to (Lx) natural(rho) = -normal(0,k*dx(rho))! Flex PDE default BCs ! start (0) line to (Lx) ! Flex PDE default BCs TIME 0 TO tmax MONITORS elevation(rho) from (0) to (Lx) PLOTS FOR t = 0 BY .01*Tmax TO tmax elevation(rho) from (0) to (Lx) summary report(VAL(rho,0)) report(VAL(rho,Lx/2)) report(VAL(rho,Lx)) END | ||
| Jared Barber (jared_barber) Member Username: jared_barber Post Number: 33 Registered: 01-2007 |
Hey, A quick note, if anything it seems you are assigning the value of rho(x,t) ((0) value(rho)) to be equal to the value of rho(x,t) (val(rho,0)) at x = 0 (and x = Lx). I think this is equivalent to saying 0 = 0 and doesn't give any information to the problem to help the solution process (I tried it on my FlexPDE with failure resulting). I'm guessing you want something more like value(rho) = 1. If I were you I would do the following: In the definitions section, add rho_init = exp(-(x)**2/(2*sigma^2))/(2*PI*sigma^2) In the initial values section change it to (just for ease in reading...doesn't really change anything in terms of solving the problem) rho = rho_init Then finally in the Boundaries section change it to: point value(rho) = value(rho_init,0) and point value(rho) = value(rho_init,Lx) Unfortunately I'm more used to 2-d probs so I don't know if you need to use point value rather than value or load rather than natural but hopefully the idea above works for you. It seemed to work for me though I could be off (i.e. I got the solutions I would expect to get for an insulated rod which is what the above bcs correspond to). Good luck, Jared TITLE 'Gaussian BC test' COORDINATES cartesian1 DEFINITIONS Lx=1e3 sigma =1 x0 = Lx/2 tmax=100 k=1 rho_init = exp(-(x)**2/(2*sigma^2))/(2*PI*sigma^2) VARIABLES rho_diri (threshold = 1e-2) rho_neu (threshold = 1e-2) INITIAL VALUES rho_diri = rho_init rho_neu = rho_init EQUATIONS rho_diri: dt(rho_diri)=k*dxx(rho_diri) rho_neu: dt(rho_neu)=k*dxx(rho_neu) BOUNDARIES region 'domain' start (0) point value(rho_diri) = val(rho_init,0) natural(rho_neu) = -normal(0,k*dx(rho_neu)) line to (Lx) point value(rho_diri) = val(rho_init,Lx) natural(rho_neu) = -normal(0,k*dx(rho_neu)) ! Dirichlet ! start (0) natural(rho) = -normal(0,k*dx(rho)) line to (Lx) natural(rho) = -normal(0,k*dx(rho))! Flex PDE default BCs ! start (0) line to (Lx) ! Flex PDE default BCs TIME 0 TO tmax MONITORS for cycle = 1 elevation(rho_diri-rho_neu) from (0) to (Lx) PLOTS FOR t = 0 BY .01*Tmax TO tmax elevation(rho_diri-rho_neu) from (0) to (Lx) summary report(VAL(rho_diri-rho_neu,0)) report(VAL(rho_diri-rho_neu,Lx/2)) report(VAL(rho_diri-rho_neu,Lx)) END | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 1120 Registered: 06-2003 |
Jared has hit the spot: value(rho)=val(rho,0) is a tautology equivalent to "rho is what it is". He correctly amended the value() conditions to POINT VALUE, which is necessary in 1D (otherwise the condition is applied throughout the line). But the Naturals need to be POINT NATURAL() as well, for the same reason. Normal() is not necessary, because in 1d there is only one direction, so you might as well use dx() instead, with proper attention to signs. If you use Normal(), it needs only a single component argument. | ||
| Fred Rodriguez (plasmon) New member Username: plasmon Post Number: 2 Registered: 12-2007 |
Thanks for your time and help Jared and Robert. With this information I am ready to move on to 2D and 3D problems. |