> > > Johan Hake wrote: >> On Wednesday 19 August 2009 19:59:00 [email protected] wrote: >>>> On Wednesday 19 August 2009 19:11:42 Garth N. Wells wrote: >>>>> Johan Hake wrote: >>>>>> On Wednesday 19 August 2009 18:50:13 Garth N. Wells wrote: >>>>>>> Johan Hake wrote: >>>>>>>> On Wednesday 19 August 2009 17:29:52 Garth N. Wells wrote: >>>>>>>>> What's going on behind the scene when I copy one vector to >>>>>>>>> another >>>>> in >>>>> >>>>>>>>> PyDOLFIN using >>>>>>>>> >>>>>>>>> u0.vector()[:] = u.vector()[:] >>>>>>>> When I add: >>>>>>>> >>>>>>>> u2 = Function(V) >>>>>>>> u2.vector()[:] = u.vector()[:] >>>>>>>> >>>>>>>> # Plot solution >>>>>>>> plot(u2) >>>>>>>> >>>>>>>> In a the Poisson demo it works fine. >>>>>>> In parallel? >>>>>> Yupp! >>>>>> >>>>>>>> This type of slice only uses the assignment operator. >>>>>>> The GenericVector assignment operator? >>>>>> Yes. This is done in __setslice__/__getslice__ in the extended >>>>>> python >>>>>> class. >>>>> I don't see then why the we end up inside >>>>> >>>>> _compare_vector_with_vector >>>>> >>>>> for the Cahn-Hilliard demo when running in parallel? >>>> Are you sure it is in >>>> >>>> _compare_vector_with_vector >>>> >>>> this should only kick in if you used '==' or some other comparison >>>> operator. >>>> >>>> When I run this demo in parallel, I get passed the assignment line but >>>> it >>>> stops when calling the solve function with the following message: >>>> >>>> Traceback (most recent call last): >>>> File "demo.py", line 86, in <module> >>>> Traceback (most recent call last): >>>> solver.solve(problem, u.vector()) >>>> RuntimeError: *** Error: MUMPS is required for parallel symbolic LU. >>>> >>>>>>> I delved into the problem and the program ends up in the function >>>>>>> _compare_vector_with_vector from dolfin_la_get_set_items.i. It uses >>>>>>> GenericVector::get and GenericVector::set. These need to be used >>>>>>> with >>>>>>> caution in parallel. >>>>>> Yes, when other type of slices are used we need to find a way to do >>>>> that >>>>> >>>>>> in parallel. I have just started digging in the parallel code, and >>>>> need >>>>> >>>>>> to look into all the changes to get familiar with it before I try to >>>>>> solve this. >>>>> This shouldn't be too hard. >>>>> >>>>> set(const double* block, uint m, const uint* rows) >>>>> >>>>> works in parallel but >>>>> >>>>> get(double* block, uint m, const uint* rows) >>>>> >>>>> doesn't yet (not too hard to fix though). The really evil functions >>>>> are >>>>> >>>>> set(const double*) >>>>> >>>>> and >>>>> >>>>> get(double*) >>>> I think I am only using >>>> >>>> set(const double* block, uint m, const uint* rows) >>>> get(double* block, uint m, const uint* rows) >>>> >>>> However GenericVector.array is using the get(double*) function, and >>>> returns >>>> some fishy numbers in the numpy.array. >>>> >>>> Johan >>> I guess anything that involves NumPy will turn into garbage in >>> parallel, >>> right? >> >> Should we issue an error, or should we try to work with the local >> values? >> > > If we use get/set(double*) very carefully and use the longer forms of > get/set more, and give an offet to numpy arrays if possible (this can be > got using the MPI::range function or GenericVector::range) could work > out ok. > > Garth >
I guess the numpy typemaps etc. assume contiguous memory (at least earlier versions of the typemaps did). Should the parallel arrays be made into such arrays ? Kent _______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
