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

[David K. <dav...@ak...>]

On Fri, Oct 20, 2017 at 12:08 PM, John Peterson <jwp...@gm...>
wrote:

>
>
> 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.
>


OK, I'll test that, and make a PR.

Question: Will the "dof_index = elem->dof_number..." line work when elem is
not owned by proc 0?

David