jothurgood / simplecfd Goto Github PK
View Code? Open in Web Editor NEWSimple (and not-so-simple) CFD solvers written in Fortran with Python plotting routines
License: MIT License
Simple (and not-so-simple) CFD solvers written in Fortran with Python plotting routines
License: MIT License
Problem with website that is meant to generate the html?
test 3 is test 2 but passed to var eta solver for eta = 1 everywhere
eta bcs are set to "fixed" but fixed eta for inhomo not implemented. I.e., I think it might be seeing eta=0 on the boundaries.
Why does it still work and not cause problems?
Current test is time-stepping a diffusion equation and doesnt give 2nd order convg.
Devise an independent test to check the order of MG itself when running in this mode is still 2nd.
Some weird debugging condition on the output has crept into the main source, clean it up
!DEBUG
if (step >= 1211) then
print *, '*** max divu before cleaning',maxval(abs(divu)*dt)
endif
e.g. the red-black relax loops in multigrid should be easy enough to put in a parallel do
as it says in title
To keep track of options used in a run
Architecture for pre-defined test problems has become a bit unwieldy.
Would like to within the main code essentially have user control / initial conditions, but a one line override that points to a specific problem module that has its own control, initial condition and diagnostic that takes precedent over the main module.
IO currently outputs simple .dat files and immediately calls a python script to plot them. This can lead to many calls to python (once for each dump). The python plotting script is substantially slower than the rest of the code.
Ideally
mainly because I like to have a terminal vertically split on a small laptop with fairly large font, so this reads well and doesn't visually break lines
In step 4 should probabally calc a time centred rho (from 0.5 * (rho^n + rho^n+1 ) to use in calc of U*
Currently uses rho^n+1 (as advect_dens is called prior, which updates the rho array to the new time level).
Should test and then (if good) implement in CCAPS, CCAPS_vardens and CCAPS_atm_const
want to test vardens on a Rayleigh Taylor instability
check this and fix if need be
make relaxation step in CCAPS more modular instead of two calls to very similar subroutine. Then the simple projection tests can be built from the same makefile in a cleaner manner, and a future multigrid solver could just be parachuted in-line into the advance_module
The need to have vface bcs is just related to applying minmod at the very edge, otherwise the two ghosts for the sln vars should be enough to directly give you everything you need for the edge vars.
Add handling of the extreme case to minmod (just revert to a simple gradient there), and then change some of the loop sizes to capture the limits appropriately.
If you run with appropriate compiler flags, should check if OOB for you.
The simplification of boundary conditions should remove a number of possible sources of error/bugs and make the code easier to use.
Write an exact R solver for the 1D Euler equations in an ideal gas
To provide exactly converged solutions against which to compare the various numerical solutions.
(later ... possibly make a split 2D extension, or even MHD if feeling brave / have time)
no longer true
in grav case, you still get some bad divergence left across the "opposite" no slip boundary to the direction of gravity (e.g. grav_y < 0 => y_max (no_slip) having this issue).
Not seen any bad behaviour but compiler debug warning
[ 50%] Building Fortran object CMakeFiles/CCAPS.dir/src/advance_module.f90.o
/Users/tflv6/Data/SimpleCFD/incompressible/CCAPS_atm_const/src/advance_module.f90:313:0:
betay_r = 0.5_num * (betacc + p0(iy+1)**(1.0_num / gamma))
Warning: 'betacc' may be used uninitialized in this function [-Wmaybe-uninitialized]
Inspecting code, seems will in fact occur when evaluating the iy=0 slice for ix>0.
This might be causing problems, although you've not spotted anything so far.
In advance_module
336 if (use_vardens) call rho_bcs ! needed for any OOB in relax and correction
but rho_bc was updated to need arguments about whether cc, x_face, y_face, etc.
Is this call just not updating the boundary?
Could be a big source of error!!!
some use di=0 for x components, d=1 for y. Others use 1 and 2. Fortran would typically use 1 and 2. Switch all to 1 and 2
mixed order of access
Add one at a time, from simple advection up, refactoring as go so all are consistent with one another
Description
Cause
the smaller errors at iy = ny-1 are a consequence of the BC choice of extrapolation of pressure through the boundary. Different BC choices / higher order interpolation / not using centred differences for evaluation of pressure gradient near the edges might improve this (I suspect not using CD but "inward difference" near the edge is the way to go so gradp is never determined based on boundary?)
the larger error at iy = ny is mainly due to rho(ix,iy=ny) not remaining =1 at the boundary (although the choice of pressure bc also effects it somewhat)
this is an unexpected result from advect_dens
Tracing it back, it emanates from the evaluation of the right (full/transversal) rho state)
705 rhoyr(ix,iy) = rhohyr(ix,iy) &
706 & - 0.5_num * dt * rho(ix,iy+1) * (macv(ix,iy+1) - macv(ix,iy)) / dy &
707 & - 0.5_num * dt / dx * ( rhohx(ix,iy+1)*macu(ix,iy+1) - &
708 & rhohx(ix-1,iy+1)*macu(ix-1,iy+1) )
Specifically the term at 706.
This is because you "expected" macv(ix,iy+1) - macv(ix,iy)
to evaluate as 0 but no slip means that macv(iy+1) = - macv(iy), so you get a non zero term
So the boundary as is doesn't really make sense with the gravity? Could either implement a zero-gradient / inoutflow boundary at the top, or implement (after this step, gradp should correctly compensate the gravity and noslip on internal flows makes more sense?)
When first wrote these was using
https://github.com/jacobwilliams/pyplot-fortran
Useful, but gives it a bit of a modular dependency and that FoBis nonsense...
In interests of portability of the whole repository, better to just write own python scripts to read *dat and call from within Fortran so that images are automatically produced during run time (if so desired)
currently not properly tested or set up for all faces and flavours.
Think its possible (0.001_num / (time+0.001_num))
should be all square rooted, and its just worked out to be fine anyway because you only simulate for small time.
!the analytical solution
do ix = 1,nx
do iy = 1,ny
xx = x(ix) - 0.5_num
yy = y(iy) - 0.5_num
r = sqrt(xx**2 + yy**2)
phi_an(ix,iy) = (phi2-phi1) * (0.001_num / (time+0.001_num)) &
& * exp(-0.25_num * r**2 / k / (0.001_num + time)) + phi1
enddo
enddo
Hard-coding fixed eta values for the MG boundary will be a hassle in the long run.
Also, current hardcoding corresponding to rho = 1 and rho = 7 for RTI isn't really correct? it should be those ^-1 . Why does the RTI mostly work? Why did you find you start to get nonconvergence issues when you correct it?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.