From: Geordie McBain <gdmcbain@fr...>  20120724 08:00:08

2012/7/23 Stephane Popinet <s.popinet@...>: > Hi Geordie, > >> What do you make of the output of this little Gerris script? I took >> out all the particles and am just accelerating the fluid. I figured >> that with an even acceleration a, provided by a GfsSource, the >> velocity should be U (t) = a t. But no, it looks like a (t + dt). > > This is clearly a (serious) bug. What happens is that the initial > velocity field is projected before the first timestep (to make sure > that it is nondivergent), this is done here in simulation_run(): > > if (sim>time.i == 0) { > gfs_approximate_projection (domain,... > > For a consistent treatment of (centered) momentum source terms, they > are included as a pressure gradient term during the projection (in > timestep.c:mac_projection()). This first projection thus includes an > acceleration term where it should not. One can check that by setting: > > Time { i = 1 end = 0.9 } > > which gets rid of the first projection. I confirm that. I hadn't known of that behaviour of the i keyword in GfsTime, but I guess it's an obvious thing to do for restarts. So in fact this issue is not at all what I thought it was, i.e. not an istart=1 issue. For comparison, the attached heating.gfs which is the analogous steady uniform GfsSource for a GfsVariableTracer produces an off by half a timestep result: %< T time: 0 first: 5.000e02 second: 5.000e02 infty: 5.000e02 T time: 0.1 first: 1.500e01 second: 1.500e01 infty: 1.500e01 T time: 0.2 first: 2.500e01 second: 2.500e01 infty: 2.500e01 T time: 0.3 first: 3.500e01 second: 3.500e01 infty: 3.500e01 T time: 0.4 first: 4.500e01 second: 4.500e01 infty: 4.500e01 T time: 0.5 first: 5.500e01 second: 5.500e01 infty: 5.500e01 T time: 0.6 first: 6.500e01 second: 6.500e01 infty: 6.500e01 T time: 0.7 first: 7.500e01 second: 7.500e01 infty: 7.500e01 T time: 0.8 first: 8.500e01 second: 8.500e01 infty: 8.500e01 T time: 0.9 first: 1.073e+154 second: 1.073e+154 infty: 1.073e+154 >% which is consistent with what you wrote about tracers in the concurrent thread. This is still a firstorder error in time though. > I can't think of a quick solution but this needs to be fixed. O.K. Thanks, Stephane. 