On Fri, Oct 20, 2017 at 10:17 AM, David Knezevic <david.kneze...@akselos.com
> wrote:
> On Fri, Oct 20, 2017 at 12:08 PM, John Peterson <jwpeter...@gmail.com>
> wrote:
>
>>
>>
>> On Fri, Oct 20, 2017 at 9:53 AM, David Knezevic <
>> david.kneze...@akselos.com> wrote:
>>
>>> On Fri, Oct 20, 2017 at 11:38 AM, John Peterson <jwpeter...@gmail.com>
>>> wrote:
>>>
>>>>
>>>>
>>>> On Fri, Oct 20, 2017 at 9:19 AM, David Knezevic <
>>>> david.kneze...@akselos.com> 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?
>
I think that's right, yes. We just do some indexing magic that doesn't
require parallel communication in DofObject::dof_number().
--
John
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
