search

clawpack / pyclaw

PyClaw is a Python-based interface to the algorithms of Clawpack and SharpClaw. It also contains the PetClaw package, which adds parallelism through PETSc.
http://www.clawpack.com/pyclaw
BSD 3-Clause "New" or "Revised" License
155 stars 96 forks source link

check qbc_backup for variable dt in solver.py #72

Closed mparsani closed 12 years ago

mparsani commented 12 years ago

I think that the way we backup the old solution (qold) with BC does not work properly. I have fixed it locally ( Not yet pushed to github.).

mandli commented 12 years ago

What exactly is the problem? Is the copy function not working in the Python?

Kyle

On Tue, Oct 4, 2011 at 4:29 AM, Matteo Parsani reply@reply.github.com wrote:

I think that the way we backup the old solution (qold) with BC does not work properly. I have fixed it locally ( Not yet pushed to github.).

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72

mparsani commented 12 years ago

Yes the copy does not work correctly. A temporary fix that I've introduced in the shallow-sphere branch when dt is variable is qold=qnew.copy('F'). This is what we also have for a constant dt. However, this fix does not work if the time is rejected because qold is overwritten. So that part should be modified. If you have already some idea let me know. I will work on it tomorrow.

Matteo

On 04/ott/2011, at 19:13, Kyle Mandlireply@reply.github.com wrote:

What exactly is the problem? Is the copy function not working in the Python?

Kyle

On Tue, Oct 4, 2011 at 4:29 AM, Matteo Parsani reply@reply.github.com wrote:

I think that the way we backup the old solution (qold) with BC does not work properly. I have fixed it locally ( Not yet pushed to github.).

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2288232

mandli commented 12 years ago

Is it this code that's the problem? 

if self.dt_variable:
    self.qbc_backup = self.qbc.copy('F')
     told = solution.t
 retake_step = False  # Reset flag

Is it then the fact we are copying the qbc array rather than the actual state q? Sorry, I am just trying to get a handle on what the problem is since this is a pretty critical piece of code.

mparsani commented 12 years ago

I am using the Strang splitting method for the source term. After the first src call with dt/2 the solution is correct (I am comparing with clawpack 4.3). Before the step2 is called qnew and qold should be the same at the beginning of a new the time step (q at the current time level) . However, this is not the case. I guess that in self.qbc_backup = self.qbc.copy('F') (which is very simple) something does not work. If you output qold and qnew inside step2, before calling flux2 you notice that.

By using qold=qnew.copy('F') I fixed it.

So, to see the problem, you can checkout the shallow-sphere branch and output qold and qnew inside step2 (just before the call to flux2).

On 4 October 2011 21:44, Kyle Mandli < reply@reply.github.com>wrote:

Is it this code that's the problem? {{{ if self.dt_variable: self.qbc_backup = self.qbc.copy('F') told = solution.t retake_step = False # Reset flag }}} Is it then the fact we are copying the qbc array rather than the actual state q? Sorry, I am just trying to get a handle on what the problem is since this is a pretty critical piece of code.

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2289772

Matteo

mandli commented 12 years ago

I may be that changing the copy at that code from a copy to a deepcopy may fix it if I understood your problem correctly.

Now that I look at this code, I am a little confused as to what's going on. The qbc array is populated with q via the state method get_qbc_from_q and then the copy is made. The algorithm looks like it only operates on the solution then but this of course is not the case. The values of qbc are the ones updated and then in homogeneous step are put directly into q instead of calling the appropriate routine from state. I don't think this is the correct behavior in any case as we should only modify the state's q array outside of the homogeneous step, more than likely in the base solver class. This actually brings up another interesting point, I don't think we need the qbc backup at all as the state.q should provide that as long as we do not modify q until we know that we have accepted the complete step.

Maybe the right way to do this is to subclass states so there is an idea of a qbc containing state and that can be passed around. The other issue is that we are doing extra copying in the non-adaptive time stepping (extra copy from qbc to q). I am going to continue to look at the code but if anyone has some rebuttal to my quick glance at the 1d routines, please feel free to better inform me.

Kyle

On Tue, Oct 4, 2011 at 1:59 PM, Matteo Parsani reply@reply.github.com wrote:

I am using the Strang splitting method for the source term. After the first src call with dt/2 the solution is correct (I am comparing with clawpack 4.3). Before the step2 is called qnew and qold should be the same at the beginning of a new the time step (q at the current time level) . However, this is not the case. I guess that in self.qbc_backup = self.qbc.copy('F') (which is very simple) something does not work. If you output qold and qnew inside step2, before calling flux2 you notice that.

By using qold=qnew.copy('F') is fix it.

So, to see the problem, you can checkout the shallow-sphere branch and output qold and qnew inside step2 (just before the call to flux2).

On 4 October 2011 21:44, Kyle Mandli < reply@reply.github.com>wrote:

Is it this code that's the problem? {{{ if self.dt_variable:    self.qbc_backup = self.qbc.copy('F')    told = solution.t retake_step = False  # Reset flag }}} Is it then the fact we are copying the qbc array rather than the actual state q?  Sorry, I am just trying to get a handle on what the problem is since this is a pretty critical piece of code.

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2289772

Matteo

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2289954

mparsani commented 12 years ago

Thanks Kyle!

I am looking at clawpack.py and solver.py and I can not understand what's wrong in the backup of qbc, i.e. in this piece of code (clawpack.py around line 592):

        if self.dt_variable:

            qold = self.qbc_backup # Solver should quarantee that 
                                               # qbc_backup will not be
                                               # changed so that it can be used in
                                               # case of step rejection.

When I look at solver. py where self.qbc_backup is defined I simply see

        if self.dt_variable:
            self.qbc_backup = self.qbc.copy('F')
            told = solution.t

This looks correct (we want to copy qbc and not simply the state q) but it does not do what one expects! If you are at time t^(n) and you want to compute the solution at t^(n+1), when you enter the function step2 for the first time (not yet rejected time steps) qold and qnew must be the same. However they aren't.

I also agree that there are other few other things that could be improved/modified.

1) state.q should be modified out of homogeneous step

2) If 1) is implemented we do not need qbc backup right.

The extra copy in the non-adaptive time stepping is there just because we want to use the same call to step2 function for bot adaptive and non-adaptive time stepping procedure. I think we could remove this extra copy but we should have different call like

        if self.dt_variable:

            q, cfl = classic2.step2(maxm,mx,my,mbc,mx,my, \
                  QOLD,QNEW,self.auxbc,dx,dy,dt,self.method,self.mthlim,cfl, \
                  self.aux1,self.aux2,self.aux3,self.work)

         else:
          q, cfl = classic2.step2(maxm,mx,my,mbc,mx,my, \
                  QNEW,QNEW,self.auxbc,dx,dy,dt,self.method,self.mthlim,cfl, \
                  self.aux1,self.aux2,self.aux3,self.work)
mandli commented 12 years ago

Have you tried a deepcopy instead of just a copy at the copy line?

Kyle

On Wed, Oct 5, 2011 at 1:22 AM, Matteo Parsani reply@reply.github.com wrote:

Thanks Kyle!

I am looking at clawpack.py and solver.py and I can not understand what's wrong in the backup of qbc, i.e. in this piece of code (clawpack.py around line 592):

           if self.dt_variable:

               qold = self.qbc_backup # Solver should quarantee that                                                   # qbc_backup will not be                                                   # changed so that it can be used in                                                   # case of step rejection.

When I look at solver. py where self.qbc_backup is defined I simply see

           if self.dt_variable:                self.qbc_backup = self.qbc.copy('F')                told = solution.t

This looks correct (we want to copy qbc and not simply the state q) but it does not do what one expects! If you are at time t^(n) and you want to compute the solution at t^(n+1), when you enter the function step2 for the first time (not yet rejected time steps) qold and qnew must be the same. However they aren't.

I also agree that there are other few other things that could be improved/modified.

1) state.q should be modified out of homogeneous step

2) If 1) is implemented we do not need qbc backup right.

The extra copy in the non-adaptive time stepping is there just because we want to use the same call to step2 function for bot adaptive and non-adaptive time stepping procedure. I think we could remove this extra copy but we should have different call like

           if self.dt_variable:

               q, cfl = classic2.step2(maxm,mx,my,mbc,mx,my, \                      QOLD,QNEW,self.auxbc,dx,dy,dt,self.method,self.mthlim,cfl, \                      self.aux1,self.aux2,self.aux3,self.work)

            else:              q, cfl = classic2.step2(maxm,mx,my,mbc,mx,my, \                      QNEW,QNEW,self.auxbc,dx,dy,dt,self.method,self.mthlim,cfl, \                      self.aux1,self.aux2,self.aux3,self.work)

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2294993

mparsani commented 12 years ago

With deepcopy the adaptive time step works fine. Thanks Kyle! Indeed, there is no need of an additional work array because the algorithm as it is now is okay. However, as pointed out, we should avoid to update state.q in clawpack.py (with state.set_q_from_qbc(mbc,self.qbc), in the homogeneous step function) and then in solver.py check whether the time step was rejected and reset state.q to qbc_backup (with state.set_q_from_qbc(self.mbc, self.qbc_backup), solver.py, line 682).

What's your opinion David?

Thanks again!

Matteo

ketch commented 12 years ago

In this discussion, it's important to understand how we use state.q and solver.qbc. The current philosophy is that state.q is the state, and qbc is just a very temporary array needed by the solver. Hence state.q should be updated whenever we finish part of a step (like the source term or the homogeneous step). This approach lends itself naturally to parallelization, but implies some small inefficiencies due to copies. This is also why the state is a parallel object but the solver is essentially not.

Kyle and Matteo, you both seem to be suggesting the following change to the code: Instead of updating q after each part of the algorithm, just populate qbc at the beginning of the timestep, never modify q during the timestep, and finally set q from qbc at the end after checking whether the step is accepted.

This is straightforward and efficient in the serial code, but won't work in the parallel code without significant changes, because the solver would need control of the DA, we'd need to change some globalToLocal() calls to localToLocal() calls, etc. I'm not saying we shouldn't do it, but you should be aware that it would require a significant refactoring.

ketch commented 12 years ago

After more consideration, I would definitely oppose the idea of making qbc a persistent object (i.e., using it throughout a step for different routines and expecting it to be up-to-date). The reason is that it makes the code much less modular. Indeed, I would really prefer to move to a multiphysics kind of model where the hyperbolic part is on an equal footing with other terms. Then some of the other physics might require different numbers of ghost cells (for instance).

mparsani commented 12 years ago

Thanks David. I see the importance of your point. I did not consider the parallelization aspect. So, personally I think that if it's not a priority we can leave the code as it is for now. The code works fine with deepcopy and now is running.

Matteo

On 07/ott/2011, at 14:24, David Ketchesonreply@reply.github.com wrote:

In this discussion, it's important to understand how we use state.q and solver.qbc. The current philosophy is that state.q is the state, and qbc is just a very temporary array needed by the solver. Hence state.q should be updated whenever we finish part of a step (like the source term or the homogeneous step). This approach lends itself naturally to parallelization, but implies some small inefficiencies due to copies. This is also why the state is a parallel object but the solver is essentially not.

Kyle and Matteo, you both seem to be suggesting the following change to the code: Instead of updating q after each part of the algorithm, just populate qbc at the beginning of the timestep, never modify q during the timestep, and finally set q from qbc at the end after checking whether the step is accepted.

This is straightforward and efficient in the serial code, but won't work in the parallel code without significant changes, because the solver would need control of the DA, we'd need to change some globalToLocal() calls to localToLocal() calls, etc. I'm not saying we shouldn't do it, but you should be aware that it would require a significant refactoring.

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2320372

mandli commented 12 years ago

I am going to move this to the pyclaw dev list as I think we have actually solved this particular bug.

Kyle

On Fri, Oct 7, 2011 at 6:41 AM, Matteo Parsani < reply@reply.github.com>wrote:

Thanks David. I see the importance of your point. I did not consider the parallelization aspect. So, personally I think that if it's not a priority we can leave the code as it is for now. The code works fine with deepcopy and now is running.

Matteo

On 07/ott/2011, at 14:24, David Ketchesonreply@reply.github.com wrote:

In this discussion, it's important to understand how we use state.q and solver.qbc. The current philosophy is that state.q is the state, and qbc is just a very temporary array needed by the solver. Hence state.q should be updated whenever we finish part of a step (like the source term or the homogeneous step). This approach lends itself naturally to parallelization, but implies some small inefficiencies due to copies. This is also why the state is a parallel object but the solver is essentially not.

Kyle and Matteo, you both seem to be suggesting the following change to the code: Instead of updating q after each part of the algorithm, just populate qbc at the beginning of the timestep, never modify q during the timestep, and finally set q from qbc at the end after checking whether the step is accepted.

This is straightforward and efficient in the serial code, but won't work in the parallel code without significant changes, because the solver would need control of the DA, we'd need to change some globalToLocal() calls to localToLocal() calls, etc. I'm not saying we shouldn't do it, but you should be aware that it would require a significant refactoring.

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2320372

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2320494

ketch commented 12 years ago

As I suspected, this issue was not related to deepcopy vs copy. I have removed the deepcopy in commit 7664c1626d94e4bd6e8 and it works fine. I will close this issue when the shallow-sphere branch is merged to master.

mandli commented 12 years ago

What ended up being the problem?

ketch commented 12 years ago

We were using the same array for two different purposes. It only mattered when using Strang splitting, which is why we didn't catch it sooner.

On Tue, Oct 11, 2011 at 9:41 PM, Kyle Mandli < reply@reply.github.com>wrote:

What ended up being the problem?

Reply to this email directly or view it on GitHub: https://github.com/clawpack/pyclaw/issues/72#issuecomment-2369483