Commit 831bfd7c authored by Alexander Danilov's avatar Alexander Danilov
Browse files

Merge branch 'master' of inmost-github

parents f3e4e583 5452e1ea
......@@ -31,7 +31,7 @@ set(HEADER inmost.h
container.hpp
io.hpp
solver_ilu2.hpp
# solver_ddpqiluc2.hpp
solver_ddpqiluc2.hpp
solver_bcgsl.hpp
solver_prototypes.hpp)
......@@ -56,7 +56,7 @@ option(USE_PARTITIONER_PARMETIS "Use ParMetis partitioner" OFF)
option(USE_PARTITIONER_ZOLTAN "Use Zoltan partitioner" OFF)
option(USE_SOLVER_PETSC "Use PETSc solver" OFF)
option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF)
option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" ON)
option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF)
option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates for automatic differentiation" OFF)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake_find/")
......
This diff is collapsed.
......@@ -1427,11 +1427,14 @@ namespace INMOST
assert(array != NULL);
array = array - beg_index;
for (IndType i = beg_index; i < end_index; ++i) new (array + i) ValType(c);
//std::cout << __FUNCTION__ << " address " << array << std::endl;
}
else array = NULL;
}
interval(const interval & other)
{
//std::cout << __FUNCTION__ << std::endl;
beg_index = other.beg_index;
end_index = other.end_index;
if( beg_index != end_index )
......@@ -1439,12 +1442,17 @@ namespace INMOST
array = static_cast<ValType *>(malloc(sizeof(ValType)*(end_index-beg_index)));
assert(array != NULL);
array = array - beg_index;
for(IndType i = beg_index; i < end_index; ++i) new (array+i) ValType(other.array[i]);
for(IndType i = beg_index; i < end_index; ++i)
{
new (array+i) ValType(other.array[i]);
}
//std::cout << __FUNCTION__ << " address " << array << std::endl;
}
else array = NULL;
}
~interval()
{
//std::cout << __FUNCTION__ << " delete address " << array << std::endl;
for(IndType i = beg_index; i < end_index; i++) (array[i]).~ValType();
if( beg_index != end_index ) free(array+beg_index);
array = NULL;
......@@ -1462,6 +1470,7 @@ namespace INMOST
assert(array != NULL);
array = array - beg_index;
for(IndType i = beg_index; i < end_index; ++i) new (array+i) ValType(other.array[i]);
//std::cout << __FUNCTION__ << " address " << array << std::endl;
}
else
{
......@@ -1518,13 +1527,18 @@ namespace INMOST
}
void set_interval_end(IndType end)
{
for(IndType i = end; i < end_index; ++i) (array[i]).~ValType();
if( end == end_index ) return;
if( beg_index != end )
{
array = static_cast<ValType *>(realloc(array+beg_index,sizeof(ValType)*(end-beg_index)));
assert(array != NULL);
array = array - beg_index;
for(IndType i = end_index; i < end; ++i) new (array+i) ValType();
ValType * array_new = static_cast<ValType *>(malloc(sizeof(ValType)*(end-beg_index)));
assert(array_new != NULL);
array_new = array_new - beg_index;
for(IndType i = beg_index; i < std::min(end,end_index); ++i) new (array_new+i) ValType(array[i]);
for(IndType i = end_index; i < end; ++i) new (array_new+i) ValType();
for(IndType i = beg_index; i < end_index; ++i) array[i].~ValType();
if( array != NULL ) free(array+beg_index);
array = array_new;
}
else
{
......@@ -1818,6 +1832,15 @@ namespace INMOST
}
}
public:
void report_addr()
{
std::cout << "stack: " << &stack << std::endl;
std::cout << "pbegin: " << pbegin << std::endl;
std::cout << "pend: " << pend << std::endl;
std::cout << "preserved: " << preserved << std::endl;
std::cout << "size: " << pend-pbegin << std::endl;
std::cout << "reserved: " << preserved-pbegin << std::endl;
}
void reserve(enumerator n)
{
//std::cout << n << std::endl;
......@@ -1890,6 +1913,7 @@ namespace INMOST
}
dynarray(const dynarray & other)
{
//std::cout << __FUNCTION__ << std::endl;
enumerator n = other.size();
preallocate(n);
for(enumerator k = 0; k < n; k++)
......
project(Discretization)
set(SOURCE discr_common.cpp
discr_tpfa.cpp
discr_tpfa2.cpp
discr_ntpfa_a.cpp
discr_ntpfa_b.cpp
discr_ntpfa_c.cpp
discr_mpfa_a.cpp
discr_mpfa_b.cpp
discr_nmpfa_a.cpp
discr_nmpfa_b.cpp
discr_mpfa_o.cpp
discr_mpfa_l.cpp
discr_mpfa_c.cpp
discr_triplet_basic.cpp
discr_basic_a.cpp
discr.h)
add_library(Discretization ${SOURCE})
if(USE_MPI)
target_link_libraries(Discretization ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(Discretization PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
\ No newline at end of file
#pragma once
#ifndef _DISCR_H
#include "../../inmost.h"
using namespace INMOST;
//TODO:
//move orphan function inside class
//introduce discr_triplet_basic class
void tensor_prod1(Storage::real * K, Storage::real v[3], Storage::real out[3]);
void tensor_prod2(Storage::real * K, Storage::real v[3], Storage::real out[3]);
void tensor_prod3(Storage::real * K, Storage::real v[3], Storage::real out[3]);
int tensor_prod(int Ktype, Storage::real_array & K,
Storage::real v[3], Storage::real out[3]);
void tensor_add1(Storage::real * K, Storage::real Kout[9]);
void tensor_add2(Storage::real * K, Storage::real Kout[9]);
void tensor_add3(Storage::real * K, Storage::real Kout[9]);
int tensor_add(int Ktype, Storage::real_array & K, Storage::real Kout[9]);
int average_tensor(int Ktype, Storage::real_array & K1,
Storage::real_array & K2, Storage::real Kout[9]);
void vec_diff(Storage::real x[3], Storage::real y[3], Storage::real out[3]);
Storage::real dot_prod(Storage::real x[3], Storage::real y[3]);
//type of the averaging point
#define AVG_NONE 0//material is homogeneous across cells - no point required
#define AVG_CENTROID 1//use centroid
#define AVG_HARMONIC 2//use harmonic averaging point
#define AVG_NONLINEAR 4//use nonlinear flux continuety to define value at the point
#define AVG_REVDIST 8//use reverse distance
#define AVG_NEUMANN 16 //use zero neumann boundary condition
#define CRITICAL_DIST_RATIO 10.0
class discr_basic
{
protected:
Automatizator * aut;
Mesh * m;
Tag K;
Tag bnd_conds;
//Tag AvgCoord, AvgCoefs, AvgStats;
int Ktype;
int error;
ElementType bnd_types;
MIDType bnd_markers, add_markers;
int rem_marker;
int have_bnd;
int allowed_types;
// returns
// 0 if point inside
// 1 if distance to point by distance to farthest node is less then CRITICAL_DIST_RATIO
// 2 otherwise
int isPointInside(Face * f, Storage::real x[3]);
//void InitializeSupportPoints(int allowed_types = AVG_NONE | AVG_CENTROID | AVG_HARMONIC | AVG_NONLINEAR | AVG_REVDIST);
int FindHarmonicPoint(Face * fKL, Cell * cK, Cell * cL,
Tag tensor, int tensor_type,
Storage::real xK[3], Storage::real xL[3],
Storage::real y[3], Storage::real & coef);
int FindBoundaryPoint(Face * fK, Cell * cK,
Tag tensor, int tensor_type,
Storage::real xK[3], Storage::real y[3]);
virtual void Build(Face * f) = 0;
public:
enum DiscrType {TPFA, MPFA, NTPFA, NMPFA};
enum SchemeType{ Linear, Nonlinear };
public:
discr_basic(Automatizator * aut, Mesh * m, std::string tensor_name);
virtual bool needBuild(Face * f) const;
virtual bool needUpdate(Face * f) const;
virtual void Init() = 0;
virtual DiscrType GetDiscrType() const = 0;
virtual SchemeType GetSchemeType() const = 0;
virtual expr Grad(const expr & param) const = 0;
virtual expr BoundaryGrad(const expr & param) const = 0;
virtual expr Interp(ElementType etype,const expr & param) const = 0;
virtual void Update() = 0;
virtual void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const; //export to ADGPRS
// by default no boundary information is provided, schemes avoid using boundary
// first marker is the marker that is used inside Automatizator to define dynamic tag domain
// if some boundary element is not used inside scheme it will not be marked by usage_marker
// second marker defines boundary faces/edges, if not set, it is computed automatically
// third parameter defines what kinds of elements may be used
void Boundary(MIDType bnd_marker = 0, ElementType bnd_types = FACE | EDGE);
// mark additional collocation points that should be introduced into the system
MIDType UnknownMarker() const { return add_markers; }
MIDType BoundaryMarker() const { return bnd_markers; }
virtual ~discr_basic();
void AllowTypes(int types) {allowed_types = types;}
};
class discr_triplet_basic : public discr_basic
{
protected:
int find_triplet(Element * r, Element * f, Face * harmonic_face, Storage::real v[3],
Storage::real dir, Storage::real area,
Element * t[3], Storage::real coefs[3], int ptypes[3],
int allowed = AVG_HARMONIC | AVG_NONLINEAR | AVG_REVDIST, bool f_necessery = false);
public:
discr_triplet_basic(Automatizator * aut, Mesh * m, std::string tensor_name) : discr_basic(aut,m,tensor_name) {}
virtual ~discr_triplet_basic();
};
class tpfa : public discr_basic
{
INMOST_DATA_ENUM_TYPE stncl;
Tag trans;
Tag elems;
void Build(Face * f);
public:
tpfa(Automatizator * aut, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return TPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); }
bool needUpdate(Face * f) const { return discr_basic::needUpdate(f); }
~tpfa();
};
class tpfa2 : public discr_basic
{
INMOST_DATA_ENUM_TYPE stncl;
Tag trans;
Tag elems;
void Build(Face * f);
public:
tpfa2(Automatizator * aut, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return TPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); }
bool needUpdate(Face * f) const { return discr_basic::needUpdate(f); }
~tpfa2();
};
class mpfa_o : public discr_basic
{
INMOST_DATA_ENUM_TYPE stncl;
Tag trans;
Tag elems;
void Build(Face * f);
public:
mpfa_o(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return MPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~mpfa_o();
};
class basic_a : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp, add;
INMOST_DATA_ENUM_TYPE trpr;
INMOST_DATA_ENUM_TYPE trpl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag trans;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
basic_a(Automatizator * out, Mesh *m, std::string tensor);
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
virtual ~basic_a();
};
class ntpfa_a : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp;
INMOST_DATA_ENUM_TYPE trpr, addr;
INMOST_DATA_ENUM_TYPE trpl, addl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag trpl_add, trpr_add;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void MarkRecursiveEdge(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
ntpfa_a(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NTPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~ntpfa_a();
};
class nmpfa_a : public basic_a
{
public:
nmpfa_a(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
~nmpfa_a();
};
class ntpfa_c : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp;
INMOST_DATA_ENUM_TYPE trpr1, trpl1, addl1, addr1;
INMOST_DATA_ENUM_TYPE trpr2, trpl2, addl2, addr2;
Tag trpl1_elems, trpl2_elems, trpr1_elems, trpr2_elems;
Tag trpl1_coefs, trpl2_coefs, trpr1_coefs, trpr2_coefs;
Tag trpl1_add, trpl2_add, trpr1_add, trpr2_add;
Tag intrp_stncl, intrp_coefs;
void MarkRecursive(Element * f);
void MarkRecursiveEdge(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
ntpfa_c(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NTPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const;
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~ntpfa_c();
};
class ntpfa_b : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp;
INMOST_DATA_ENUM_TYPE trpr, addr;
INMOST_DATA_ENUM_TYPE trpl, addl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag trpl_add, trpr_add;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void MarkRecursiveEdge(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
ntpfa_b(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NTPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~ntpfa_b();
};
class nmpfa : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp, add;
INMOST_DATA_ENUM_TYPE trpr;
INMOST_DATA_ENUM_TYPE trpl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag addval;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
nmpfa(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Interp(FACE,param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~nmpfa();
};
class mpfa_t : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE scheme, interp;
Tag stncl, coefs, interp_stncl, interp_coefs;
void MarkRecursive(Element * f);
void Build(Face * f);
public:
mpfa_t(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~mpfa_t();
};
class mpfa_y : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE scheme, interp;
Tag stncl, coefs, interp_stncl, interp_coefs;
void MarkRecursive(Element * f);
void Build(Face * f);
public:
mpfa_y(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~mpfa_y();
};
#endif
\ No newline at end of file
#include "discr.h"
#include <iomanip>
#define EPS 1e-12
typedef Storage::real vec[3];
static Storage::real lengthS(vec x) { return dot_prod(x, x); }
static Storage::real length(vec x) { return sqrt(lengthS(x)); }
static Storage::real normalize(Storage::real x[3]) { Storage::real r = length(x); if (r > 0.0) x[0] /= r, x[1] /= r, x[2] /= r; return r; }
basic_a::basic_a(Automatizator * aut, Mesh * m, std::string tensor)
:discr_triplet_basic(aut, m, tensor)
{
trpl_elems = m->CreateTag("A_TRP_L", DATA_REFERENCE, FACE, NONE, 4);
trpr_elems = m->CreateTag("A_TRP_R", DATA_REFERENCE, FACE, NONE, 4);
trpl_coefs = m->CreateTag("A_COEF_L", DATA_REAL, FACE, NONE, 4);
trpr_coefs = m->CreateTag("A_COEF_R", DATA_REAL, FACE, NONE, 4);
trans = m->CreateTag("A_TRANS", DATA_REAL, FACE, NONE, 1);
trpr = aut->RegisterStencil("a_right_triplet", trpr_elems, trpr_coefs);
trpl = aut->RegisterStencil("a_left_triplet", trpl_elems, trpl_coefs);
add = aut->RegisterStaticTag(trans);
intrp_stncl = m->CreateTag("INTRP_STNCL", DATA_REFERENCE, FACE | EDGE | NODE, NONE);
intrp_coefs = m->CreateTag("INTRP_COEFS", DATA_REAL, FACE | EDGE | NODE, NONE);
dwn = aut->RegisterStencil("back_cell", get_l);
upw = aut->RegisterStencil("front_cell", get_r);
intrp = aut->RegisterStencil("interp",intrp_stncl,intrp_coefs,add_markers);
}
void basic_a::MarkRecursive(Element * f)
{
if( f->GetElementType() == FACE && f->GetMarker(bnd_markers) )
{
Storage::reference_array lrtrp[2];
lrtrp[0] = f->ReferenceArrayDF(trpl_elems);
lrtrp[1] = f->ReferenceArrayDF(trpr_elems);
for(int i = 0; i < 2; i++)
for(Storage::reference_array::iterator it = lrtrp[i].begin(); it != lrtrp[i].end(); ++it) if( (*it) != NULL )
{
if( (*it)->GetElementType() != CELL && !(*it)->GetMarker(add_markers))
{
(*it)->SetMarker(add_markers);
(*it)->SetMarker(mark_next);
}
}
for(int i = 0; i < 2; i++)
for(Storage::reference_array::iterator it = lrtrp[i].begin(); it != lrtrp[i].end(); ++it) if( (*it) != NULL )
{
if( (*it)->GetMarker(mark_next) )
{
(*it)->RemMarker(mark_next);