3

i have been struggling to compute a particular instance of cylindrical 3D heat equation.

enter image description here Here is my code : 

(*parameter*)
\[Alpha] = 0.1;
h = 1;


eq = D[u[t, r, \[Theta], z], t] - 
   D[u[t, r, \[Theta], z], r, r] + (1/r) D[u[t, r, \[Theta], z], r] + 
   1/(\[Alpha]^2 r^2) D[u[t, r, \[Theta], z], \[Theta], \[Theta]] + 
   D[u[t, r, \[Theta], z], z, z];

(*initial and boundary conditions*)

ic = {u[0, r, \[Theta], z] == Exp[-r]};
bc = {u[t, 10^-6, \[Theta], z] == 1, u[t, 1, \[Theta], z] == 0, 
   PeriodicBoundaryCondition[u[t, r, \[Theta], z], \[Theta] == 0, 
    TranslationTransform[{0, 2*Pi*\[Alpha], 0}]], 
   PeriodicBoundaryCondition[u[t, r, \[Theta], z], z == 0, 
    TranslationTransform[{0, 0, h}]]};


(*solution*)

sol = NDSolve[{eq == 
     NeumannValue[0, r == 10^-6] + NeumannValue[1, r == 1], ic, bc}, 
   u[t, r, \[Theta], z], {t, 0, 10}, {r, 10^-6, 1}, {\[Theta], 0, 
    2*Pi}, {z, 0, h}];



(*plot*)

Manipulate[
  Plot3D[sol, {r, 0, 1}, {\[Theta], 0, 2 Pi}], {t, 0, 10}, {z, 0, 1}];

I have corrected the error regarding the definition of eq. I also corrected the error regarding my periodic boundary conditions. I previously made a mistake

The period regarding the \[Theta] variable should have been 2*Pi*\[Alpha] instead of just 2*Pi. And the period regarding the z variable would just be the height of my cylinder which i called h.

edit : So it seems that the solution would be to apply TranslationTransform to only the 3 variables of space and not the time space vector. It gets rid of the mapping error but now I must face another one

NDSolve::femnodpbc: DirichletCondition can not be present on the target boundary of a PeriodicBoundaryConditon.

It appears that I must examine where the case where both the pbc and regular boundary conditions are both conflicting.

ConfuzzledStudent
  • 105
  • 9
  • 2
    The message means that 1/\[Alpha]\.b2r\.b2 does not have a value defined. Something is odd with your typesetting. Try: 1/(\[Alpha]^2 r^2) – user21 Jun 01 '20 at 10:56
  • Also note that Derivative[0, 1, 0, 0][u][t, 1, \[Theta], z] is not going to wok. Use NeumannValue for that. Or, alternatively, get rid of the PeriodicBoundaryConditon and make use of the u[t, r, 0, z]==u[t, r, 2 Pi, z] syntax. The also replace the second PBC in the same manner. – user21 Jun 01 '20 at 11:00
  • i understand, i will edit it right now. Can you explain why either replacing the derivatives by NeumanValues OR getting rid of the periodicboundaryconditions will work ? – ConfuzzledStudent Jun 01 '20 at 11:02
  • 2
    See the section What Triggers the Use of the Finite Element Method. Any objects mentioned there will trigger the FEM (which is fine) but then you need to think about The Relation between NeumannValue and Boundary Derivatives. The other alternative is to make use of the tensor product grid and for that you need to remove the FEM triggers. – user21 Jun 01 '20 at 11:14
  • i edited my code as you suggested – ConfuzzledStudent Jun 01 '20 at 11:38
  • 2
    Where did you see the usage of NeumannValue[Derivative...]? – user21 Jun 01 '20 at 11:42
  • i guess my original question was answered, i'm not sure if i should close this post and make a new one since i have a new problem :/ – ConfuzzledStudent Jun 01 '20 at 16:11
  • @ConfuzzledStudent Do you need really to solve this problem? – Alex Trounev Jun 01 '20 at 16:20
  • @AlexTrounev yes, this is very important for my research and i have been wrestling with a much bigger (more parameters) diffusion equation than this one right here, for 3 months now. Solving this simplified version of the heat diffusion equation i'm currently working on will allow me to understand better how to deal with the bigger equation, and therefore will finally allow me to progress with my research. – ConfuzzledStudent Jun 01 '20 at 16:27
  • @ConfuzzledStudent Is this bigger equation linear or nonlinear? As a rule numerical algorithm not transferable from one problem to another. So, it could be better you show your real problem. – Alex Trounev Jun 01 '20 at 16:39
  • @AlexTrounev it's the same initial and boundary conditions, what is different is an additional derivative D[u[t, r, \[Theta], z], z, \[Theta]] in the equation. Unfortunately i wish i could tell you more but i have been asked not to reveal any more than the simplified equation. I apologize for not being able to tell you more. – ConfuzzledStudent Jun 01 '20 at 16:48
  • 2
    @ConfuzzledStudent Ok! Then what is the reason you use 3 periodic boundary conditions in one direction? Is it a typo or you try to get very specific solution? – Alex Trounev Jun 01 '20 at 19:02
  • my first boundary condition is for z : u[t, r, \[Theta], 0] == u[t, r, \[Theta], 2 Pi*\[Alpha]]

    whereas the second one is for the variable theta : u[t, r, 0, z] == u[t, r, 2 Pi, z].

    There should only be two periodic boundary conditions. Did i perhaps make a third one ?

     – ConfuzzledStudent Jun 01 '20 at 19:42
  • …Why do you modify mol in that way? You're making almost the same mistake in every edit. In short: If you want to use FiniteElement method, then you cannot have Derivative in your b.c.s, you need to re-express them with NeumannValue. Please calm down, read the answers you found carefully, and stop coding blindly. – xzczd Jun 02 '20 at 10:22
  • @xzczd like i said i'm a total beginner when it comes to mathematica. I added the part about mol and its use in NDSolve because i saw that it was used in an answer to another question regarding non numerical derivative values at 0. And sure enough it got rid of the problem i had at t=0 (or 10^-6 actually). Truth be told i do not really know what it does. – ConfuzzledStudent Jun 02 '20 at 11:10
  • @xzczd Also, i got the part that if i wanted to use FEM, I needed to replace derivatives with NeumannValue or alternatively get rid of the pbc commands and write stuff like u[t, r, 0, z]==u[t, r, 2 Pi, z] instead, as user21 pointed.

    But when i did, i encountered the error NDSolve::fembcdepderiv: Derivatives of dependent variables in boundary conditions are not supported with the Finite Element Method in this version of NDSolve.. which i have no idea how to solve.

     – ConfuzzledStudent Jun 02 '20 at 11:11
  • As said above, you're coding blindly. The usage of NeumannValue isn't like that. Please check the document of it carefully by pressing F1. "Truth be told i do not really know what it does. " Then as asked above, why do you modify the mol function that way? I've never defined a mol function involving FiniteElement. Do notice mol is a helper function for the TensorProduceGrid method. – xzczd Jun 02 '20 at 11:25
  • @xzczd i understand. I will remove it. When you say that the usage of NeumannValue is wong, are you perhaps referring to my second edit ? Is so i will try my best to edit it correctly. – ConfuzzledStudent Jun 02 '20 at 11:34
  • 1
    Yes, the NeumannValue[Derivative...] is just wrong. – xzczd Jun 02 '20 at 11:45
  • @AlexTrounev yes i see where i added a third pbc, ty for pointing that out – ConfuzzledStudent Jun 02 '20 at 11:48
  • @xzczd i corrected the mistake on NeumannValue (hopefully), deleted a typo for the pbc. Btw could you tell me if the TranslationTransform is used correctly ?

    Otherwise i'm dealing with a very peculiar issue right now (see edit) but i'll try to understand the answer i found first.

     – ConfuzzledStudent Jun 02 '20 at 12:11
  • 1
  • If one of the b.c. is Derivative[0, 1, 0, 0][u][t, 1, θ, z] == 1, then the NeumannValue is still incorrect. 2. You're setting both PeriodicBoundaryCondition in θ direction, which is undoubted wrong. Once again, please calm down, and check carefully.
    • – xzczd Jun 02 '20 at 12:20
  • Thank you for pointing out the pbc error, i realise now that's also what @AlexTrounev pointed out and i missed it completely.

    Now that the pbc error is fixed, it leaves us with the NeumannValue. Would this change work correctly ? My bc was indeed Derivative[0, 1, 0, 0][u][t, 1, θ, z] == 1

    here is the edit : sol = NDSolve[{eq == NeumannValue[0, r == 10^-6] + NeumannValue[1, r == 1], ic, bc}, u[t, r, \[Theta], z], {t, 0, 10}, {r, 10^-6, 1}, {\[Theta], 0, 2*Pi}, {z, 0, 1}];

     – ConfuzzledStudent Jun 02 '20 at 13:00
  • 1
  • How did you define eq? 2. The question has become rather messy now, please consider removing some of the unsuccessful trial to keep it clean. 3.It's better to show us the problem in traditional math notation so we can check easier.
    • – xzczd Jun 02 '20 at 21:14
  • 1
    When defining eq you've already added ==, so you've now writen down an equation in the form (… = …) = NeumannValue[…], which is clearly wrong. If you still don't understand what I mean, then once again, please check the document of NeumannValue and see what's the correct way to add it to the equation. – xzczd Jun 03 '20 at 01:45
  • 2
    As @xzczd suggested writing down the equation seems a good idea. Here is a code in the correct syntax for the tensor product grid method: ic = {u[0, r, \[Theta], z] == Exp[-r]}; bc = {u[t, 10^-6, \[Theta], z] == 1, u[t, 1, \[Theta], z] == 0, u[t, r, 0, z] == u[t, r, 2 Pi, z], u[t, r, \[Theta], 0] == u[t, r, \[Theta], 2*Pi], Derivative[0, 1, 0, 0][u][t, 10^-6, \[Theta], z] == 0, Derivative[0, 1, 0, 0][u][t, 1, \[Theta], z] == 1 }; sol = NDSolve[{eq == 0, ic, bc}, u[t, r, \[Theta], z], {t, 0, 10}, {r, 10^-6, 1}, {\[Theta], 0, 2*Pi}, {z, 0, 1}]; – user21 Jun 03 '20 at 05:58
  • 2
    but that still has issues. If you want to make use of FEM you'd need to write the equation in inactive form. I suggest you try to get the tensor product grid (TPG) method to work first – user21 Jun 03 '20 at 05:59
  • @xzczd i see what you mean there, thank you for pointing that out. – ConfuzzledStudent Jun 03 '20 at 10:08
  • Alright, i will try to write down the equation as well as the boundary conditions as clearly as possible. – ConfuzzledStudent Jun 03 '20 at 10:14
  • I wrote down the equation and the initial/boundary conditions so that it may be as clear as possible. I also corrected some previous mistake i made regarding the role played by the variable \[Alpha] . – ConfuzzledStudent Jun 03 '20 at 13:42
  • @user21 I assume FEM would require me to write my equation into a formal PDE judging by the page that comes up when i click on the three dot next to the error. I did not really understand what it means, but i would rather not my equation too much.

    I guess i should look at how to use the TPG method then.

     – ConfuzzledStudent Jun 03 '20 at 13:48
  • I do not see where you have specified a z boundary condition. You specify it to be periodic, but with no value specified at z=0 or z=h so it could be anything. – Bill Watts Jun 03 '20 at 23:15
  • Also your bc's show no dependence on $\theta$ other than periodic, so your solution will be constant in $\theta$. This might confuse NDSolve if you leave $\theta$ in your pde with no pertinent bc's so that term should be removed. If your more complicated equation requires $\theta$ dependence and you want to include it in this problem then you should specify some $\theta$ dependence in either your initial condition or at your r = 1 boundary. – Bill Watts Jun 04 '20 at 00:30
  • @Bill The topic about periodic b.c. is a bit involved. In short, though in many materials it's not explicitly mentioned, usually derivative in $z$ direction is also required to be periodic, but currently PeriodicBoundaryCondition doesn't behave like this, for more info see the discussion in this post. Currently the simplest way to get the traditional periodic b.c. (in regular domain of course) is to use the old good TensorProductGrid. Then, 5th b.c. looks strange, why is the period $2\pi \alpha$? – xzczd Jun 04 '20 at 02:58
  • @xzczd Well my point was from just a differential equation point of view rather than a Mathematica point of view, if there is no $\theta$ dependence in the boundary conditions, there will be no $\theta$ dependence in the solution. – Bill Watts Jun 04 '20 at 05:44
  • @Bill I'm commenting on your first comment i.e. b.c. in $z$ direction. My point is a traditional periodic b.c. in $z$ direction is enough, we don't need additional b.c. at $z=0$ or $z=h$. – xzczd Jun 04 '20 at 05:57
  • @xzczd Regarding the 5th condition : u(t,r,θ,z) =u(t,r,θ+ 2πα) The parameter α represents an angular deficit (or surplus if α>1). So, the rotation around z, in this case, cannot reach up to 2π.....Sorry that's the best i can explain with my english skills. – ConfuzzledStudent Jun 04 '20 at 10:07
  • Then how can the b.c. in $\theta$ direction be periodic? – xzczd Jun 04 '20 at 10:27
  • There shouldn't be any problem as long as α is a rational number right ? – ConfuzzledStudent Jun 04 '20 at 13:01
  • So this cylinder does not go all the way around? Sounds like it has a wedge of size $2\pi-2\pi \alpha$. – Bill Watts Jun 04 '20 at 18:50
  • @BillWatts That is indeed what alpha represents. Compared to the normal cylindrical diffusion equation, you remove a wedge off of the cylinder and try to link the discontinuous ends/faces....Not sure if that's clear. – ConfuzzledStudent Jun 04 '20 at 19:02
  • The mapping error is no more, but another took its place... – ConfuzzledStudent Jun 10 '20 at 10:05