Skip to content

Commit

Permalink
Added arguments_dict and plugged it into AddedMass
Browse files Browse the repository at this point in the history
  • Loading branch information
JohanMabille authored and cekees committed Apr 24, 2020
1 parent ff9eabe commit 270e257
Show file tree
Hide file tree
Showing 6 changed files with 308 additions and 41 deletions.
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,8 @@ matlab_setup.done
dist/
tmp/

timeHistory.txt

# Generated files

proteus/BoundaryConditions.cpp
Expand Down
48 changes: 44 additions & 4 deletions proteus/mprans/AddedMass.h
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <iostream>
#include "CompKernel.h"
#include "ModelFactory.h"
#include "ArgumentsDict.h"

#include "xtensor-python/pyarray.hpp"

Expand All @@ -15,7 +16,8 @@ namespace proteus
{
public:
virtual ~cppAddedMass_base(){}
virtual void calculateResidual(//element
virtual void calculateResidual(arguments_dict& args) = 0;
/*virtual void calculateResidual(//element
xt::pyarray<double>& mesh_trial_ref,
xt::pyarray<double>& mesh_grad_trial_ref,
xt::pyarray<double>& mesh_dof,
Expand Down Expand Up @@ -52,7 +54,7 @@ namespace proteus
xt::pyarray<double>& Aij,
int added_mass_i,
xt::pyarray<double>& barycenters,
xt::pyarray<int>& flags_rigidbody)=0;
xt::pyarray<int>& flags_rigidbody)=0;*/
virtual void calculateJacobian(//element
xt::pyarray<double>& mesh_trial_ref,
xt::pyarray<double>& mesh_grad_trial_ref,
Expand Down Expand Up @@ -225,7 +227,7 @@ namespace proteus
}
}

void calculateResidual(//element
/*void calculateResidual(//element
xt::pyarray<double>& mesh_trial_ref,
xt::pyarray<double>& mesh_grad_trial_ref,
xt::pyarray<double>& mesh_dof,
Expand Down Expand Up @@ -262,8 +264,46 @@ namespace proteus
xt::pyarray<double>& Aij,
int added_mass_i,
xt::pyarray<double>& barycenters,
xt::pyarray<int>& flags_rigidbody)
xt::pyarray<int>& flags_rigidbody)*/
void calculateResidual(arguments_dict& args)
{
xt::pyarray<double>& mesh_trial_ref = args.m_darray["mesh_trial_ref"];
xt::pyarray<double>& mesh_grad_trial_ref = args.m_darray["mesh_grad_trial_ref"];
xt::pyarray<double>& mesh_dof = args.m_darray["mesh_dof"];
xt::pyarray<int>& mesh_l2g = args.m_iarray["mesh_l2g"];
xt::pyarray<double>& dV_ref = args.m_darray["dV_ref"];
xt::pyarray<double>& u_trial_ref = args.m_darray["u_trial_ref"];
xt::pyarray<double>& u_grad_trial_ref = args.m_darray["u_grad_trial_ref"];
xt::pyarray<double>& u_test_ref = args.m_darray["u_test_ref"];
xt::pyarray<double>& u_grad_test_ref = args.m_darray["u_grad_test_ref"];
//element boundary
xt::pyarray<double>& mesh_trial_trace_ref = args.m_darray["mesh_trial_trace_ref"];
xt::pyarray<double>& mesh_grad_trial_trace_ref = args.m_darray["mesh_grad_trial_trace_ref"];
xt::pyarray<double>& dS_ref = args.m_darray["dS_ref"];
xt::pyarray<double>& u_trial_trace_ref = args.m_darray["u_trial_trace_ref"];
xt::pyarray<double>& u_grad_trial_trace_ref = args.m_darray["u_grad_trial_trace_ref"];
xt::pyarray<double>& u_test_trace_ref = args.m_darray["u_test_trace_ref"];
xt::pyarray<double>& u_grad_test_trace_ref = args.m_darray["u_grad_test_trace_ref"];
xt::pyarray<double>& normal_ref = args.m_darray["normal_ref"];
xt::pyarray<double>& boundaryJac_ref = args.m_darray["boundaryJac_ref"];
//physics
int nElements_global = args.m_iscalar["nElements_global"];
int nElementBoundaries_owned = args.m_iscalar["nElementBoundaries_owned"];
xt::pyarray<int>& u_l2g = args.m_iarray["u_l2g"];
xt::pyarray<double>& u_dof = args.m_darray["u_dof"];
xt::pyarray<double>& q_rho = args.m_darray["q_rho"];
int offset_u = args.m_iscalar["offset_u"];
int stride_u = args.m_iscalar["stride_u"];
xt::pyarray<double>& globalResidual = args.m_darray["globalResidual"];
int nExteriorElementBoundaries_global = args.m_iscalar["nExteriorElementBoundaries_global"];
xt::pyarray<int>& exteriorElementBoundariesArray = args.m_iarray["exteriorElementBoundariesArray"];
xt::pyarray<int>& elementBoundaryElementsArray = args.m_iarray["elementBoundaryElementsArray"];
xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.m_iarray["elementBoundaryLocalElementBoundariesArray"];
xt::pyarray<int>& elementBoundaryMaterialTypesArray = args.m_iarray["elementBoundaryMaterialTypesArray"];
xt::pyarray<double>& Aij = args.m_darray["Aij"];
int added_mass_i = args.m_iscalar["added_mass_i"];
xt::pyarray<double>& barycenters = args.m_darray["barycenters"];
xt::pyarray<int>& flags_rigidbody = args.m_iarray["flags_rigidbody"];
for(int eN=0;eN<nElements_global;eN++)
{
for (int k=0;k<nQuadraturePoints_element;k++)
Expand Down
75 changes: 39 additions & 36 deletions proteus/mprans/AddedMass.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
from proteus.TransportCoefficients import TC_base
from proteus.SubgridError import SGE_base
from proteus.ShockCapturing import ShockCapturing_base
from . import cAddedMass
from . import cAddedMass, cArgumentsDict
from mpi4py import MPI

class NumericalFlux(proteus.NumericalFlux.ConstantAdvection_Diffusion_SIPG_exterior):
Expand Down Expand Up @@ -668,43 +668,46 @@ def getResidual(self, u, r):
u[self.offset[0] + self.stride[0] * dofN] = g(self.dirichletConditionsForceDOF[0].DOFBoundaryPointDict[dofN], self.timeIntegration.t)
self.setUnknowns(u)
self.Aij[:,:,self.added_mass_i]=0.0
self.addedMass.calculateResidual( # element
self.u[0].femSpace.elementMaps.psi,
self.u[0].femSpace.elementMaps.grad_psi,
self.mesh.nodeArray,
self.mesh.elementNodesArray,
self.elementQuadratureWeights[('u', 0)],
self.u[0].femSpace.psi,
self.u[0].femSpace.grad_psi,
self.u[0].femSpace.psi,
self.u[0].femSpace.grad_psi,

argDict = cArgumentsDict.ArgumentsDict()
argDict.darray['mesh_trial_ref'] = self.u[0].femSpace.elementMaps.psi
argDict.darray['mesh_grad_trial_ref'] = self.u[0].femSpace.elementMaps.grad_psi
argDict.darray['mesh_dof'] = self.mesh.nodeArray
argDict.iarray['mesh_l2g'] = self.mesh.elementNodesArray
argDict.darray['dV_ref'] = self.elementQuadratureWeights[('u', 0)]
argDict.darray['u_trial_ref'] = self.u[0].femSpace.psi
argDict.darray['u_grad_trial_ref'] = self.u[0].femSpace.grad_psi
argDict.darray['u_test_ref'] = self.u[0].femSpace.psi
argDict.darray['u_grad_test_ref'] = self.u[0].femSpace.grad_psi
# element boundary
self.u[0].femSpace.elementMaps.psi_trace,
self.u[0].femSpace.elementMaps.grad_psi_trace,
self.elementBoundaryQuadratureWeights[('u', 0)],
self.u[0].femSpace.psi_trace,
self.u[0].femSpace.grad_psi_trace,
self.u[0].femSpace.psi_trace,
self.u[0].femSpace.grad_psi_trace,
self.u[0].femSpace.elementMaps.boundaryNormals,
self.u[0].femSpace.elementMaps.boundaryJacobians,
argDict.darray['mesh_trial_trace_ref'] = self.u[0].femSpace.elementMaps.psi_trace
argDict.darray['mesh_grad_trial_trace_ref'] = self.u[0].femSpace.elementMaps.grad_psi_trace
argDict.darray['dS_ref'] = self.elementBoundaryQuadratureWeights[('u', 0)]
argDict.darray['u_trial_trace_ref'] = self.u[0].femSpace.psi_trace
argDict.darray['u_grad_trial_trace_ref'] = self.u[0].femSpace.grad_psi_trace
argDict.darray['u_test_trace_ref'] = self.u[0].femSpace.psi_trace
argDict.darray['u_grad_test_trace_ref'] = self.u[0].femSpace.grad_psi_trace
argDict.darray['normal_ref'] = self.u[0].femSpace.elementMaps.boundaryNormals
argDict.darray['boundaryJac_ref'] = self.u[0].femSpace.elementMaps.boundaryJacobians
# physics
self.mesh.nElements_global,
self.mesh.nElementBoundaries_owned,
self.u[0].femSpace.dofMap.l2g,
self.u[0].dof,
self.coefficients.q_rho,
self.offset[0], self.stride[0],
r,
self.mesh.nExteriorElementBoundaries_global,
self.mesh.exteriorElementBoundariesArray,
self.mesh.elementBoundaryElementsArray,
self.mesh.elementBoundaryLocalElementBoundariesArray,
self.mesh.elementBoundaryMaterialTypes,
self.Aij,
self.added_mass_i,
self.barycenters,
self.flags_rigidbody)
argDict.iscalar['nElements_global'] = self.mesh.nElements_global
argDict.iscalar['nElementBoundaries_owned'] = self.mesh.nElementBoundaries_owned
argDict.iarray['u_l2g'] = self.u[0].femSpace.dofMap.l2g
argDict.darray['u_dof'] = self.u[0].dof
argDict.darray['q_rho'] = self.coefficients.q_rho
argDict.iscalar['offset_u'] = self.offset[0]
argDict.iscalar['stride_u'] = self.stride[0]
argDict.darray['globalResidual'] = r
argDict.iscalar['nExteriorElementBoundaries_global'] = self.mesh.nExteriorElementBoundaries_global
argDict.iarray['exteriorElementBoundariesArray'] = self.mesh.exteriorElementBoundariesArray
argDict.iarray['elementBoundaryElementsArray'] = self.mesh.elementBoundaryElementsArray
argDict.iarray['elementBoundaryLocalElementBoundariesArray'] = self.mesh.elementBoundaryLocalElementBoundariesArray
argDict.iarray['elementBoundaryMaterialTypesArray'] = self.mesh.elementBoundaryMaterialTypes
argDict.darray['Aij'] = self.Aij
argDict.iscalar['added_mass_i'] = self.added_mass_i
argDict.darray['barycenters'] = self.barycenters
argDict.iarray['flags_rigidbody'] = self.flags_rigidbody
self.addedMass.calculateResidual(argDict)
# sum of residual should be zero
# but it is sometimes not exactly zero with certain meshes in parallel
# hack
Expand Down
122 changes: 122 additions & 0 deletions proteus/mprans/ArgumentsDict.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,122 @@
#include "pybind11/pybind11.h"
#include "pybind11/stl_bind.h"

#define FORCE_IMPORT_ARRAY
#include "ArgumentsDict.h"

namespace xt
{
#if defined(__GNUC__) && !defined(__clang__)
namespace workaround
{
inline void complex_allocator()
{
std::allocator<int> ai;
std::allocator<double> ad;
}
}
#endif
}

namespace py = pybind11;

namespace proteus
{
template <class K, class T>
void bind_pyarray_dict(py::class_<pyarray_dict<K, T>>& cl, const std::string& name)
{
using map_type = pyarray_dict<K, T>;
using key_type = typename map_type::key_type;
using mapped_type = typename map_type::mapped_type;

cl.def(py::init<>());

cl.def("__setitem__",
[](map_type& m, const key_type& k, mapped_type& v)
{
auto it = m.find(k);
if(it != m.end())
{
it->second = std::move(v);
}
else
{
m.emplace(k, std::move(v));
}
}
);

cl.def("__repr__",
[name](map_type& m)
{
std::ostringstream oss;
oss << name << '{';
bool f = false;
for(const auto& kv: m)
{
if(f)
{
oss << ", ";
}
oss << kv.first <<": " << kv.second;
f = true;
}
oss << '}';
return oss.str();
},
"Return the canonical representation of this map."
);

cl.def("__bool__",
[](map_type& m) -> bool
{
return !m.empty();
},
"Check wether the map is nonempty"
);

cl.def("__getitem__",
[](map_type& m, const key_type& k) -> mapped_type&
{
auto it = m.find(k);
if(it == m.end())
{
throw std::runtime_error(detail::key_error_message(k));
}
return it->second;
}
);
}
}

PYBIND11_MAKE_OPAQUE(proteus::scalar_dict<double>);
PYBIND11_MAKE_OPAQUE(proteus::scalar_dict<int>);

PYBIND11_MODULE(cArgumentsDict, m)
{
using proteus::pyarray_dict;
using proteus::scalar_dict;
using proteus::arguments_dict;

xt::import_numpy();

using dpyarray_dict = pyarray_dict<std::string, double>;
using ipyarray_dict = pyarray_dict<std::string, int>;

py::class_<dpyarray_dict> dad(m, "DArrayDict");
proteus::bind_pyarray_dict(dad, "DArrayDict");

py::class_<ipyarray_dict> iad(m, "IArrayDict");
proteus::bind_pyarray_dict(iad, "IArrayDict");

py::bind_map<scalar_dict<double>>(m, "DScalarDict");
py::bind_map<scalar_dict<int>>(m, "IScalarDict");

py::class_<arguments_dict>(m, "ArgumentsDict")
.def(py::init<>())
.def_readwrite("darray", &arguments_dict::m_darray)
.def_readwrite("iarray", &arguments_dict::m_iarray)
.def_readwrite("dscalar", &arguments_dict::m_dscalar)
.def_readwrite("iscalar", &arguments_dict::m_iscalar);
}

Loading

0 comments on commit 270e257

Please sign in to comment.