Re: [Libmesh-users] How to use ExodusII_IO::copy_elemental_solution in parallel?

[John P. <jwp...@gm...>]

On Fri, Oct 20, 2017 at 9:53 AM, David Knezevic <dav...@ak...>
wrote:

> On Fri, Oct 20, 2017 at 11:38 AM, John Peterson <jwp...@gm...>
> wrote:
>
>>
>>
>> On Fri, Oct 20, 2017 at 9:19 AM, David Knezevic <
>> dav...@ak...> wrote:
>>
>>> I'm using ExodusII_IO::copy_elemental_solution to copy element data to
>>> a new System object, and this works fine in serial using code like the
>>> below:
>>>
>>> ----------------------------------------------
>>>
>>>   Mesh mesh(comm);
>>>   ExodusII_IO exo_io(mesh);
>>>
>>>   if (comm.rank() == 0)
>>>   {
>>>     exo_io.read(path_to_mesh_file);
>>>   }
>>>   MeshCommunication().broadcast(mesh);
>>>
>>>   mesh.prepare_for_use();
>>>
>>>   EquationSystems es(mesh);
>>>   ExplicitSystem& materials_system =
>>>     es->add_system<ExplicitSystem> ("data");
>>>   materials_system.add_variable ("data", CONSTANT, MONOMIAL);
>>>   es->init();
>>>
>>>   exo_io.copy_elemental_solution(es, "data", "data");
>>>
>>> ----------------------------------------------
>>>
>>>
>>> However, in parallel I hit the error "ERROR, ExodusII file must be
>>> opened for reading before copying an elemental solution!" at the start of
>>> copy_elemental_solution.
>>>
>>> As far as I know (e.g. based on namebased_io.C) we should only call
>>> "read" on proc 0, so I don't see how to make this work in parallel, since
>>> in that case exo_io is only "opened for reading" on proc 0.
>>>
>>> Does anyone have any suggestions on this? (John, I gather that you were
>>> the author of most of the ExodusII reader stuff, so I thought you might
>>> know the answer?)
>>>
>>
>> Hmm... I think copy_elemental_solution() is only designed to work in
>> serial, so could we just guard the parts the parts of that function that
>> read from the file in "if (rank==0)"?
>>
>> In addition, we would only error out if the file is not open for reading
>> on proc 0...
>>
>
> OK, that makes sense. I think we could just wrap the whole method (up to
> solution->close) in "if (rank==0)", and remove the "if" around
> system.solution->set so that all values are set on proc 0 and then
> communicated when we close the vector. Do you agree?
>

Yeah, that inner if-test around solution->set doesn't need to be there, and
solution->close() is the first collective I see.

-- 
John