Hello, I hope you are well. I am trying to solve the heat equation in parabolic cylindrical coordinates, using the finite difference method for the spatial part of the equation and Runge-Kutta 4 for the temporal part. But when running the code most of my temperature values give me zero and some other random values and I don’t know what I’m doing wrong. A first inspection tells me that perhaps my variables are not well declared, as well as my subroutines and therefore the program fails, but overall I don’t know what is failing in my code. I come here asking for help, since I am not that skilled in Fortran and I think I could receive some help here.
Here’s my code. I share the code on github because at the moment I don’t know how to share the code directly here
Hello. You mention cylindrical co-ordinates, but your notation (dx, dy) does not reflect this as one would expect something like (dr, dz) instead. Nevermind - my suggestion is to start by making tfirst the simpler case using cartesian co-ordinates in 2D or even 1D, maybe also usng just Euler integration instead of RK. Once you have the working code you can expand it to the cylindrical coordinates and RK integration.
It seems to me that the formula on line 139 does not make sense from the mathematical point of view: To calculate the temporal derivative of the temperatureI would expect something like:
“tempout(i2,j2)=a*(tempxx/dx/dx + tempyy/dy/dy)”. The values of actual position x(i2) and y(j2) do not seem to be needed here.
In addition to @grofz’s comments, I do not see (in a quick glance) what your boundary conditions are .In the solve routine the values on x=1 and y=ny are overwritten with the values at x=2 (zero because of the starting condition) and y=ny-1 (ditto).
We may also comment on the source code itself - two quick comments:
You only need to specify implicit none in the module, not for every subroutine in that module.
Use a consistent indentation. That makes it much easier to recognise the structure (do-loops, if-blocks, …)
No, I’m not using cylindrical coordinates, but parabolic cylindrical coordinates, that are related to cartesian coordinates in the following way
X = 1/2(u^2 + v^2)
Y = u*v
Z = Z
And the heat equation looks like this
and that’s why the line 139 have that equation, but I’m just solving it in 2D. I’m using X and Y variables for simplicity.
Unless I’m missing something, in your RK4 routine, you appear to be passing the variables x1, x2, x3 to solve before they are defined.
Shouldn’t it be something like:
x1=tempo+ktemp1*dt2
call solve(t+dt2,x1,tempout)!Paso 2 de 4
ktemp2=tempout
!deallocate(tempout)
!deallocate(x1)
x2=tempo+ktemp2*dt2
call solve(t+dt2,x2,tempout) !Paso 3 de 4
ktemp3=tempout
!deallocate(tempout)
!deallocate(x2)
x3=tempo+ktemp3*dt
call solve(t+dt,x3,tempout) !Paso 4 de 4
ktemp4=tempout
Remember, Fortran does not automatically initialize variables to zero. You have to make sure every variable that appears on the RHS of an assignment has a legitimate value before you try to use it or you will get NaNs etc.
- two quick comments: