Skip to content

Rework coupling model #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jgiret wants to merge 3 commits into
base: reworkCouplingModel
Choose a base branch
from
Open

Rework coupling model #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3b07fa7
Create test for GEMS toy problem
Oct 5, 2020
82b4f67
Add test for disciplines, chain and MDA
Nov 30, 2020
d5f3c06
Modification on coupling cpp and hpp files
Nov 30, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions python-wrapper/build.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
python setup.py build_ext --inplace --bitpit-path=/home/jc.giret/bitpit/ --madeleine-path=/home/jc.giret/SOURCES/madeleine_parallel/build/src/ --extensions-source=coupling.pyx --metis-path=/home/jc.giret/.local/lib/ --lapack-path=/home/jc.giret/.conda/envs/madeleine_py38/lib/ --mpi-include-path=/home/jc.giret/.conda/envs/madeleine_py38/lib/ --petsc-path=/home/jc.giret/.conda/envs/madeleine_py38/lib/
2 changes: 2 additions & 0 deletions python-wrapper/env.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
export LD_LIBRARY_PATH=/home/jc.giret/bitpit/lib/:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH=/home/jc.giret/SOURCES/madeleine_parallel/build/src/:$LD_LIBRARY_PATH
135 changes: 84 additions & 51 deletions python-wrapper/gems_wrapper.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,23 +7,39 @@
# initial documentation
# :author: Francois Gallard
# OTHER AUTHORS - MACROSCOPIC CHANGES
from __future__ import print_function
from numpy import ones
import numpy as np
from petsc4py import PETSc

from gems.core.discipline import MDODiscipline
from gems.parallel.api import create_execution_context, create_user_partition
from gems.parallel.core.mpi_manager import get_world_comm
from gems_mpi_plugins.api import create_execution_context, create_user_partition
from gems_mpi_plugins.core.mpi_manager import MPIManager
import coupling

PETSC_DETERMINE = PETSc.DETERMINE

class ToySphereDiscipline(MDODiscipline):

def __init__(self, name, inputs, outputs, mesh_file,
neutral_mesh_file, sphere_radius=1.0, sphere_neutral_radius=1.0,
sphere_thickness=0.001, is_inner_sphere=True,
source_intensity=1.0, source_direction=None,
n_cpus=1, kernel=1):
class ToySphereDiscipline(MDODiscipline):
def __init__(
self,
name,
inputs,
outputs,
mesh_file,
neutral_mesh_file,
sphere_radius=1.0,
sphere_neutral_radius=1.0,
sphere_thickness=0.001,
is_inner_sphere=True,
source_intensity=10.0,
source_direction=None,
thermalDiffusivityCoefficient=0.5,
emissivity=0.5,
infinityTemperature=400.0,
n_cpus=1,
kernel=1,
):

if source_direction is None:
source_direction = [1.0, 0.0, 0.0]
Expand All @@ -37,7 +53,7 @@ def __init__(self, name, inputs, outputs, mesh_file,

super(ToySphereDiscipline, self).__init__(name)

self.input_grammar.initialize_from_data_names([inputs[0], 'r'])
self.input_grammar.initialize_from_data_names([inputs[0], "r"])
self.output_grammar.initialize_from_data_names([outputs[0]])

# Creation of the execution context and register
Expand All @@ -49,55 +65,63 @@ def __init__(self, name, inputs, outputs, mesh_file,
local_size = 0
all_sizes = None
if self.execution_context.is_alive():
cell_indices_per_rank = coupling.Py_computeMeshFilePartitioning(neutral_mesh_file,
comm)
cell_indices_per_rank = coupling.Py_computeMeshFilePartitioning(
neutral_mesh_file.encode(), comm
)
local_indices = [i for i in cell_indices_per_rank if i == comm.rank]
local_size = len(local_indices)
all_sizes = comm.allgather(local_size)

# Broadcast all_sizes to all the ranks of comm_world
comm_world = get_world_comm()
comm_world = MPIManager().main_comm
root = self.execution_context.comm_world_ranks[0]
all_sizes = comm_world.bcast(all_sizes, root=root)

obj_variables = [0] * n_cpus
obj_variables[0] = 1

variables_sizes = {inputs[0]: all_sizes,
outputs[0]: all_sizes,
"r": obj_variables}
variables_sizes = {
inputs[0]: all_sizes,
outputs[0]: all_sizes,
"r": obj_variables,
}

partition = create_user_partition(self.execution_context,
variables_sizes)
partition = create_user_partition(self.execution_context, variables_sizes)

# Init the discipline only if the ranks belong to the discipline
# subcommunicator
partition.initialize_default_inputs(1.)
partition.initialize_default_inputs(1.0)

if self.execution_context.is_alive():
if comm.rank == 0:
self.default_inputs['r'] = ones(1)*sphere_radius
self.default_inputs["r"] = ones(1) * sphere_radius

if self.execution_context.is_alive():
self.mesh_coupling = coupling.Py_MeshCoupling(inputs,
outputs,
name,
comm)

self.mesh_coupling.initialize(mesh_file,
neutral_mesh_file,
sphere_radius,
sphere_neutral_radius,
sphere_thickness,
is_inner_sphere,
source_intensity,
source_direction,
cell_indices_per_rank,
self.kernel)
inputs_encode = [inp.encode() for inp in inputs]
outputs_encode = [out.encode() for out in outputs]
self.mesh_coupling = coupling.Py_MeshCoupling(
inputs_encode, outputs_encode, name.encode(), comm
)

self.mesh_coupling.initialize(
mesh_file.encode(),
neutral_mesh_file.encode(),
sphere_radius,
sphere_neutral_radius,
sphere_thickness,
is_inner_sphere,
source_intensity,
source_direction,
thermalDiffusivityCoefficient,
emissivity,
infinityTemperature,
cell_indices_per_rank,
self.kernel,
)

def _run(self):
input_vector = self.local_data[self.inputs[0]]
r = self.local_data['r']
r = self.local_data["r"]

comm = self.execution_context.comm
r = comm.bcast(r, root=0)
Expand All @@ -111,8 +135,9 @@ def close(self):
self.mesh_coupling.close()

def _compute_jacobian(self, inputs=None, outputs=None):
inputs = self.get_input_data_names()
outputs = self.get_output_data_names()

inputs = list(self.get_input_data_names())
outputs = list(self.get_output_data_names())

# Compute global number of neutral cell
nofNeutralLocalCells = self.mesh_coupling.getNeutralMeshSize()
Expand All @@ -122,38 +147,46 @@ def _compute_jacobian(self, inputs=None, outputs=None):
# Initialize Jacobian matrices
self._init_jacobian(with_zeros=True)
self.jac[outputs[0]] = {}
# mat = PETSc.Mat().createAIJ(nofNeutralGlobalCells,nofNeutralGlobalCells,comm=comm) #Non-preallocated, malloc at each setVelue
matInputs = PETSc.Mat().create(comm=comm)
matInputs.setSizes(((nofNeutralLocalCells, PETSC_DETERMINE), (nofNeutralLocalCells, PETSC_DETERMINE)))
matInputs.setSizes(
(
(nofNeutralLocalCells, PETSC_DETERMINE),
(nofNeutralLocalCells, PETSC_DETERMINE),
)
)
matInputs.setType("dense")
matInputs.setUp()

for i in range(nofNeutralLocalCells):
rowId,colIds,values = self.mesh_coupling.extractOutputInputJacobianRow(i)
rowId, colIds, values = self.mesh_coupling.extractOutputInputJacobianRow(i)
rowIds = [rowId]
#CAVEAT: petsc4py accepts only int32 by default. bitpit indices are long integers. Cast is possible but very large meshes are not feasible
matInputs.setValues(rowIds,colIds,values, addv=1)
# CAVEAT: petsc4py accepts only int32 by default. bitpit indices are long integers. Cast is possible but very large meshes are not feasible
matInputs.setValues(rowIds, colIds, values, addv=1)

matInputs.assemblyBegin()
matInputs.assemblyEnd()
self.jac[outpus[0]][inputs[0]] = matInputs
viewer = PETSc.Viewer().createASCII("Jac.dat",mode=1,comm=comm)
viewer.view(self.jac[outpus[0]][inputs[0]])
self.jac[outputs[0]][inputs[0]] = matInputs
viewer = PETSc.Viewer().createASCII("CouplingJac.dat", mode=1, comm=comm)
viewer.view(self.jac[outputs[0]][inputs[0]])

matControl = PETSc.Mat().create(comm=comm)
matControl.setSizes(((nofNeutralLocalCells, PETSC_DETERMINE), (1, 1)))
matControl.setType("dense")
matControl.setUp()

for i in range(nofNeutralLocalCells):
rowId,colIds,values = self.mesh_coupling.extractOutputControlJacobianRow(i)
rowId, colIds, values = self.mesh_coupling.extractOutputControlJacobianRow(
i
)
rowIds = [rowId]
#CAVEAT: petsc4py accepts only int32 by default. bitpit indices are long integers. Cast is possible but very large meshes are not feasible
mat.setValues(rowIds,colIds,values, addv=1)
colIds = [1]
# CAVEAT: petsc4py accepts only int32 by default. bitpit indices are long integers. Cast is possible but very large meshes are not feasible
print("Jac r", i, values)
matControl.setValues(rowIds, colIds, values, addv=1)

matControl.assemblyBegin()
matControl.assemblyEnd()
self.jac[outputs[0]]['r'] = matControl
self.jac[outputs[0]]["r"] = matControl

viewer = PETSc.Viewer().createASCII("RadiusJac.dat",mode=1,comm=comm)
viewer.view(self.jac[outputs[0]]['r'])
viewer = PETSc.Viewer().createASCII("RadiusJac.dat", mode=1, comm=comm)
viewer.view(self.jac[outputs[0]]["r"])
105 changes: 77 additions & 28 deletions python-wrapper/test_chain.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,46 +7,95 @@
# initial documentation
# :author: Francois Gallard
# OTHER AUTHORS - MACROSCOPIC CHANGES
from __future__ import print_function, unicode_literals
import sys, petsc4py

petsc4py.init(sys.argv)
import unittest
from numpy import ones
from math import pi
from numpy import ones, full, mean
from numpy.random import random

from gems.parallel.core.mpi_manager import MPIManager, get_world_comm
from gems.parallel.core.parallel_chain import MDOParallelMPIChain
from gems_mpi_plugins.core.mpi_manager import MPIManager
from gems_mpi_plugins.core.parallel_chain import MDOParallelMPIChain
from gems_wrapper import ToySphereDiscipline

COMM = get_world_comm()
COMM = MPIManager().main_comm
SIZE = MPIManager().main_comm.size

class TestGEMSWrapper(unittest.TestCase):

def test_basic(self):
MPIManager().clear(2)
def test_basic():
MPIManager().clear(1)

mesh_file = "../examples/data/unitSphere5.stl"
neutral_mesh_file = "../examples/data/unitSphere4.stl"
mesh_file = "../examples/data/unitSphere5.stl"
neutral_mesh_file = "../examples/data/unitSphere4.stl"

toy1 = ToySphereDiscipline("Sphere1", ["Forces"], ["Pressure"], mesh_file,
neutral_mesh_file, sphere_radius=1.0,
n_cpus=2)
toy_inner = ToySphereDiscipline(
"Sphere1",
["T_out"],
["T_in"],
mesh_file,
neutral_mesh_file,
sphere_radius=1.0,
n_cpus=1,
is_inner_sphere=True,
source_intensity=1.0,
source_direction=None,
thermalDiffusivityCoefficient=0.0,
emissivity=1e-6,
infinityTemperature=350.0,
)

toy2 = ToySphereDiscipline("Sphere1", ["Pressure"], ["Forces"], mesh_file,
neutral_mesh_file, sphere_radius=1.0,
n_cpus=2)
# neutral_mesh_size = toy_inner.mesh_coupling.getNeutralMeshSize()

chain = MDOParallelMPIChain([toy1, toy2])
default_inputs = chain.default_inputs
inputs = None
if chain.execution_context.is_rank_on_mpi_group():
r = ones(1)
size_forces = default_inputs['Forces'].shape[0]
forces = ones(size_forces)*10.
inputs = {'r': r, 'Forces': forces}
toy_outer = ToySphereDiscipline(
"Sphere1",
["T_in"],
["T_out"],
mesh_file,
neutral_mesh_file,
sphere_radius=1.0,
n_cpus=1,
is_inner_sphere=False,
source_intensity=30.0,
source_direction=None,
thermalDiffusivityCoefficient=0.25,
emissivity=0.1,
infinityTemperature=400.0,
)

out = chain.execute(inputs)
mesh_file = "../examples/data/unitSphere5.stl"
neutral_mesh_file = "../examples/data/unitSphere4.stl"

chain = MDOParallelMPIChain([toy_outer, toy_inner])
default_inputs = chain.default_inputs
inputs = None

if chain.execution_context.is_rank_on_mpi_group():
r = full(1, 10.0)
neutral_size = default_inputs["T_in"].shape[0]
t_array = full(neutral_size, 300.0) + random(neutral_size) * 10.0
t_out = {"T_out": t_array}
inputs = {"r": r}
inputs.update(t_out)

out = chain.execute(inputs)

for i in range(50):
out = chain.execute(out)

if chain.execution_context.is_rank_on_mpi_group():
target = 500.0
obj = 0.5 * mean((out["T_in"] - target) ** 2)
obj2 = 0.5 * mean((out["T_in"] - target) ** 2) * 4 * pi * r[0] ** 2
print("Obj is", obj, obj2)

# COMM.Barrier()
# if toy_inner.execution_context.is_rank_on_mpi_group():
# toy_inner.close()
# if toy_outer.execution_context.is_rank_on_mpi_group():
# toy_outer.close()

if chain.execution_context.is_rank_on_mpi_group():
print out

if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
test_basic()
Loading