Problem restarting karman.c example from dump file

[Louis Steytler]

Hi All,

When restarting from a dump file in the karman example (http://basilisk.fr/src/examples/karman.c), I get the following error:

/home/steytler/basilisk/src/embed-tree.h:280: refine_embed_linear: Assertion `((double *) ((((Tree *)grid)->L[point.level-1]->m[(point.i+2)/2+child.x] [(point.j+2)/2+0]) + sizeof(Cell)))[cs.i] && ((double *) ((((Tree *)grid)->L[point.level-1]->m[(point.i+2)/2+0] [(point.j+2)/2+child.y]) + sizeof(Cell)))[cs.i]' failed.

The original karman.c file does not include restarting from dump files or writing of dump files, so I modified the init event and included an event to dump files as shown at the bottom.

Any ideas why this error happens when restarting from a dump file? (Starting the computation from scratch produces no errors.)

In another problem, the starting cylinder test (http://basilisk.fr/src/test/starting.c), dump files are created and used to restart, but in this case does not produce any errors.

Could this be a bug or am I missing something?

Any advice would be much appreciated.

Thanks,

Louis

------------------------------------------

event init (t = 0)
{

  /**
  The domain is the intersection of a channel of width unity and a
  circle of diameter 0.125. */
 if (!restore ("restart")) {

    vertex scalar phi[];
    foreach_vertex() {
    phi[] = intersection (0.5 - y, 0.5 + y);
    phi[] = intersection (phi[], sq(x) + sq(y) - sq(0.125/2.));
    }
    boundary ({phi});
    fractions (phi, cs, fs);

  }
  else {

      vertex scalar phi[];
      foreach_vertex() {
        phi[] = intersection (0.5 - y, 0.5 + y);
        phi[] = intersection (phi[], sq(x) + sq(y) - sq(0.125/2.));
      }
      boundary ({phi});
      fractions (phi, cs, fs);

  }
  /**
  We set the initial velocity field. */
 
  foreach()
    u.x[] = cs[] ? 1. : 0.;
}

event snapshots (i += 10)
{
  char name [80];
  sprintf (name, "dump-%d", i);
  dump (name);
}

[Antoon van Hooft]

Hallo Louis,

Could you share the dump file as well?

Antoon

[Louis Steytler]

Hi Antoon,

Attached are my restart and C file.

Thanks very much,

Louis

[Louis Steytler]

Hi Antoon,

When I change the line

phi[] = intersection (0.50 - y, 0.50 + y);

to

phi[] = intersection (0.505 - y, 0.505 + y);

I get this error when restarting:

steytle1@steytle1:~/BASILISK/examples-Aug-2019/karman$ ./karman
21 0.101329 0 0
22 0.103498 2 1
23 0.106214 2 1
24 0.109442 2 1
*** Error in `./karman': munmap_chunk(): invalid pointer: 0x000000000106ab00 ***
Aborted (core dumped)
avconv version 9.20-6:9.20-0ubuntu0.14.04.2+fdkaac, Copyright (c) 2000-2014 the Libav developers
  built on Feb 13 2017 12:16:45 with gcc 4.8 (Ubuntu 4.8.4-2ubuntu1~14.04.3)
[ppm @ 0x1362920] [IMGUTILS @ 0x7ffddd8cff20] Picture size 1x0 is invalid
[image2pipe @ 0x13583c0] Estimating duration from bitrate, this may be inaccurate
Input #0, image2pipe, from 'pipe:':
  Duration: N/A, bitrate: N/A
    Stream #0.0: Video: ppm, rgb24, 1x0, 25 tbn
Output #0, mpegts, to 'vort.mp4':
Output file #0 does not contain any stream

However, when I remove the movie event it works.

Any ideas why this is?

Thanks again,

Louis

[Antoon van Hooft]

Hallo,

I believe the issue is similar as here:

http://basilisk.fr/sandbox/Antoonvh/embed_and_refine.c

For your first problem, you can work-around the issue using the following `init` event:

event init (t = 0) {

  for (scalar s in {f, u, g})
    s.prolongation = refine_injection;
  restore ("restart");
  for (scalar s in {f, u, g})
    s.prolongation = refine_embed_linear;

 
  vertex scalar phi[];
  foreach_vertex() {
    phi[] = intersection (0.5 - y, 0.5 + y);
    phi[] = intersection (phi[], sq(x) + sq(y) - sq(0.125/2.));
  }
  boundary ({phi});
  fractions (phi, cs, fs);

  boundary (all);
}

 

Antoon

[Stephane Popinet]

At the moment "dump" does not save face fields. In particular it does
not save the "fs" embedded volume fraction field. So this field needs to
be reinitialised after a "restore".

At some point, I will modify dump so that all fields are dumped. But
this is not trivial.

cheers,

Stephane

[Louis Steytler]

Antoon, Stephane,

That works, thanks very much!

Louis

--
You received this message because you are subscribed to the Google Groups "basilisk-fr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to basilisk-fr...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/ab33eab6-555a-9d7f-c2f1-62075acb6777%40basilisk.fr.

[Alessandro Della Pia]

Dear Louis,

I am working on "karman.c" test and I found this conversation on the group.

Would you mind to explain me the physical meaning of the following code lines?

u.n[embed] = fabs(y) > 0.25 ? neumann(0.) : dirichlet(0.); /* 1 */
u.t[embed] = fabs(y) > 0.25 ? neumann(0.) : dirichlet(0.);  /* 2 */

I mean, if I have to impose no-slip boundary condition on cylinder surface,

maybe I can just write:

u.n[embed] = dirichlet(0.); /* 3 */
u.t[embed] = dirichlet(0.); /* 4 */

Running the simulation replacing  lines 1-2 with lines  3-4 , results seem to be identical.

What I am not understanding about embedded boundary conditions?

Thanks in advance for your attention.

Best Regards,

Alessandro Della Pia

[Alessandro Della Pia]

up

[Louis Steytler]

Hi Alessandro,

Lines 1-2 should enforce a no-slip (Dirichlet) condition on the cylinder surface, and a free-slip (neumann) on the side walls.

Lines 3-4  should enforce no-slip (Dirichlet) at both the cylinder surface and on the side walls.

To explain why there seems to be no difference between the two cases, see this discussion:

https://groups.google.com/forum/?utm_medium=email&utm_source=footer#!searchin/basilisk-fr/louis%7Csort:date/basilisk-fr/SLnKb8U0BXs/JR9HrhvzEwAJ

Hope that helps,

Louis

--
You received this message because you are subscribed to the Google Groups "basilisk-fr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to basilisk-fr...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/504ae427-43fb-4d13-a239-131952ff29e8%40googlegroups.com.

[Alessandro Della Pia]

Hi Louis,

thanks a lot for your reply. The link that you have sent me has been very useful o understand the question.

I have just a last question:

why

u.n[embed] = fabs(y) > 0.25 ? neumann(0.) : dirichlet(0.);

u.t[embed] = fabs(y) > 0.25 ? neumann(0.) : dirichlet(0.);

and not

u.n[embed] = fabs(y) > 0.0625 ? neumann(0.) : dirichlet(0.);
u.t[embed] = fabs(y) > 0.0625 ? neumann(0.) : dirichlet(0.);

?

I mean, cylinder diameter is 0.125 = 2 * 0.0625 and cylinder centre is at y = 0.

Am I misunderstandig something about embedded boundary conditions?

Thanks in advance,

Alessandro.

To unsubscribe from this group and stop receiving emails from it, send an email to basil...@googlegroups.com.

[Louis Steytler]

The side walls are beyond 0.25.

To unsubscribe from this group and stop receiving emails from it, send an email to basilisk-fr...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/99aed6d4-adcf-4ed8-82fc-aeb8aba4dd36%40googlegroups.com.

[Alessandro Della Pia]

So, considering that the side walls are at y = + 0.5 and y = -0.5 ,

I could also write for example fabs(y) > 0.3, or fabs(y) > 0.4, obtaining the

same results? 

To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/99aed6d4-adcf-4ed8-82fc-aeb8aba4dd36%40googlegroups.com.