Commit 5452e1ea by Kirill Terekhov

Version update

some changes
parent 6c87951b
......@@ -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/")
......
......@@ -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)
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
#include "discr.h"
#include "discr.h"
#define EPS 1e-8
void mpfa_o::Build(Face * f)
{
if (needBuild(f))
{
}
}
mpfa_o::mpfa_o(Automatizator * out, Mesh *m, std::string tensor) :discr_basic(out,m, tensor)
{
trans = m->CreateTag("MPFA_TRANS_COEF", DATA_REAL, FACE, NONE, 2);
elems = m->CreateTag("MPFA_TRANS_ELEM", DATA_REFERENCE, FACE, NONE, 2);
stncl = aut->RegisterStencil("mpfa", elems, trans);
}
void mpfa_o::Init()
{
if (!have_bnd) Boundary(0, FACE);
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
{
adjacent<Face> faces = it->getFaces();
Storage::real * mat = new Storage::real[faces.size()*faces.size()];
int k = 0;
for(adjacent<Face>::iterator f = faces.begin(); f != faces.end(); f++)
{
if( f->GetMarker(bnd_markers) )
{
}
else
{
}
k++;
}
}
}
expr mpfa_o::Grad(const expr & param) const
{
return stencil(stncl, param);
}
expr mpfa_o::Interp(ElementType etype,const expr & param) const
{
return 0;
}
void mpfa_o::Export(std::ostream & fdout, Storage::real trans_scale, Storage::real vol_scale) const
{
discr_basic::Export(fdout,trans_scale,vol_scale);
Tag unknown_id = m->CreateTag("ID",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
Tag row = m->CreateTag("ROW",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
int cnt = 0, cntsupp = 0;
INMOST_DATA_ENUM_TYPE idnum = 0;
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
it->IntegerDF(unknown_id) = idnum++;
if( bnd_conds.isValid() )
for(Mesh::iteratorElement it = m->BeginElement(FACE); it != m->EndElement(); ++it)
if(it->HaveData(bnd_conds)) it->IntegerDF(unknown_id) = idnum+cntsupp++;
cnt = 0;
fdout << "MPFACONNS" << std::endl;
for (Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
{
Cell * r0 = f->BackCell(), *r1 = f->FrontCell();
if(needBuild(&*f))
{
if( f->GetMarker(bnd_markers) && (!bnd_conds.isValid() || !f->HaveData(bnd_conds))) continue;
fdout << f->ReferenceArrayDF(elems)[0]->IntegerDF(unknown_id) << " ";
fdout << f->ReferenceArrayDF(elems)[1]->IntegerDF(unknown_id) << " ";
fdout << f->RealArrayDF(trans)[0]*trans_scale << std::endl;
f->IntegerDF(row) = cnt++;
}
}
fdout << "/" << std::endl;
if( cntsupp )
{
fdout << "SUPPCONNS" << std::endl;
fdout << cntsupp << std::endl;
for (Mesh::iteratorElement f = m->BeginElement(FACE); f != m->EndElement(); ++f)
{
if( f->GetMarker(bnd_markers) && f->HaveData(bnd_conds) )
{
fdout << "--bnd " << ElementTypeName(f->GetElementType()) << "# " << f->LocalID() << std::endl;
fdout << 2 << " "; //boundary type
fdout << f->IntegerDF(row) << " "; //corresponding entry in connection list
// \alpha C + \beta dC/dn = \gamma
if( bnd_conds.isValid() && f->HaveData(bnd_conds) )
{
Storage::real_array bnd = f->RealArray(bnd_conds);
fdout << bnd[0] << " "; // \alpha
fdout << bnd[1] << " "; // \beta
fdout << bnd[2] << " "; // \gamma
}
else
{
fdout << 0.0 << " "; // \alpha
fdout << 1.0 << " "; // \beta
fdout << 0.0 << " "; // \gamma
}
fdout << std::endl;
}
}
fdout << "/" << std::endl;
}
}
void mpfa_o::Update()
{
for (INMOST_DATA_INTEGER_TYPE id = 0; id < m->MaxLocalID(FACE); ++id)
{
Face * f = m->FaceByLocalID(id);
if (f != NULL)
{
if (needBuild(f))
{
Cell * r0 = f->BackCell();
Cell * r1 = f->FrontCell();
if (r0->New() || r1->New()) Build(f);
}
}
}
}
mpfa_o::~mpfa_o()
{
m->DeleteTag(elems);
m->DeleteTag(trans);
}
\ No newline at end of file
#include "discr.h"
#include "discr.h"
#define EPS 1e-8
void mpfa_o::Build(Face * f)
{
if (needBuild(f))
{
}
}
mpfa_o::mpfa_o(Automatizator * out, Mesh *m, std::string tensor) :discr_basic(out,m, tensor)
{
trans = m->CreateTag("MPFA_TRANS_COEF", DATA_REAL, FACE, NONE, 2);
elems = m->CreateTag("MPFA_TRANS_ELEM", DATA_REFERENCE, FACE, NONE, 2);
stncl = aut->RegisterStencil("mpfa", elems, trans);
}
void mpfa_o::Init()
{
if (!have_bnd) Boundary(0, FACE);
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
{
adjacent<Face> faces = it->getFaces();
Storage::real * mat = new Storage::real[faces.size()*faces.size()];
int k = 0;
for(adjacent<Face>::iterator f = faces.begin(); f != faces.end(); f++)
{
if( f->GetMarker(bnd_markers) )
{
}
else
{
}
k++;
}
}
}
expr mpfa_o::Grad(const expr & param) const
{
return stencil(stncl, param);
}
expr mpfa_o::Interp(ElementType etype,const expr & param) const
{
return 0;
}
void mpfa_o::Export(std::ostream & fdout, Storage::real trans_scale, Storage::real vol_scale) const
{
discr_basic::Export(fdout,trans_scale,vol_scale);
Tag unknown_id = m->CreateTag("ID",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
Tag row = m->CreateTag("ROW",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
int cnt = 0, cntsupp = 0;
INMOST_DATA_ENUM_TYPE idnum = 0;
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
it->IntegerDF(unknown_id) = idnum++;
if( bnd_conds.isValid() )
for(Mesh::iteratorElement it = m->BeginElement(FACE); it != m->EndElement(); ++it)
if(it->HaveData(bnd_conds)) it->IntegerDF(unknown_id) = idnum+cntsupp++;
cnt = 0;
fdout << "MPFACONNS" << std::endl;
for (Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
{
Cell * r0 = f->BackCell(), *r1 = f->FrontCell();
if(needBuild(&*f))
{
if( f->GetMarker(bnd_markers) && (!bnd_conds.isValid() || !f->HaveData(bnd_conds))) continue;
fdout << f->ReferenceArrayDF(elems)[0]->IntegerDF(unknown_id) << " ";
fdout << f->ReferenceArrayDF(elems)[1]->IntegerDF(unknown_id) << " ";
fdout << f->RealArrayDF(trans)[0]*trans_scale << std::endl;
f->IntegerDF(row) = cnt++;
}
}
fdout << "/" << std::endl;
if( cntsupp )
{
fdout << "SUPPCONNS" << std::endl;
fdout << cntsupp << std::endl;
for (Mesh::iteratorElement f = m->BeginElement(FACE); f != m->EndElement(); ++f)
{
if( f->GetMarker(bnd_markers) && f->HaveData(bnd_conds) )
{
fdout << "--bnd " << ElementTypeName(f->GetElementType()) << "# " << f->LocalID() << std::endl;
fdout << 2 << " "; //boundary type
fdout << f->IntegerDF(row) << " "; //corresponding entry in connection list
// \alpha C + \beta dC/dn = \gamma
if( bnd_conds.isValid() && f->HaveData(bnd_conds) )
{
Storage::real_array bnd = f->RealArray(bnd_conds);
fdout << bnd[0] << " "; // \alpha
fdout << bnd[1] << " "; // \beta
fdout << bnd[2] << " "; // \gamma
}
else
{
fdout << 0.0 << " "; // \alpha
fdout << 1.0 << " "; // \beta
fdout << 0.0 << " "; // \gamma
}
fdout << std::endl;
}
}
fdout << "/" << std::endl;
}
}
void mpfa_o::Update()
{
for (INMOST_DATA_INTEGER_TYPE id = 0; id < m->MaxLocalID(FACE); ++id)
{
Face * f = m->FaceByLocalID(id);
if (f != NULL)
{
if (needBuild(f))
{
Cell * r0 = f->BackCell();
Cell * r1 = f->FrontCell();
if (r0->New() || r1->New()) Build(f);
}
}
}
}
mpfa_o::~mpfa_o()
{
m->DeleteTag(elems);
m->DeleteTag(trans);
}
\ No newline at end of file
#include "discr.h"
#include "discr.h"
#define EPS 1e-8
void mpfa_o::Build(Face * f)
{
if (needBuild(f))
{
}
}
mpfa_o::mpfa_o(Automatizator * out, Mesh *m, std::string tensor) :discr_basic(out,m, tensor)
{
trans = m->CreateTag("MPFA_TRANS_COEF", DATA_REAL, FACE, NONE, 2);
elems = m->CreateTag("MPFA_TRANS_ELEM", DATA_REFERENCE, FACE, NONE, 2);
stncl = aut->RegisterStencil("mpfa", elems, trans);
}
void mpfa_o::Init()
{
if (!have_bnd) Boundary(0, FACE);
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
{
adjacent<Face> faces = it->getFaces();
Storage::real * mat = new Storage::real[faces.size()*faces.size()];
int k = 0;
for(adjacent<Face>::iterator f = faces.begin(); f != faces.end(); f++)
{
if( f->GetMarker(bnd_markers) )
{
}
else
{
}
k++;
}
}
}
expr mpfa_o::Grad(const expr & param) const
{
return stencil(stncl, param);
}
expr mpfa_o::Interp(ElementType etype,const expr & param) const
{
return 0;
}
void mpfa_o::Export(std::ostream & fdout, Storage::real trans_scale, Storage::real vol_scale) const
{
discr_basic::Export(fdout,trans_scale,vol_scale);
Tag unknown_id = m->CreateTag("ID",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
Tag row = m->CreateTag("ROW",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
int cnt = 0, cntsupp = 0;
INMOST_DATA_ENUM_TYPE idnum = 0;
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
it->IntegerDF(unknown_id) = idnum++;
if( bnd_conds.isValid() )
for(Mesh::iteratorElement it = m->BeginElement(FACE); it != m->EndElement(); ++it)
if(it->HaveData(bnd_conds)) it->IntegerDF(unknown_id) = idnum+cntsupp++;
cnt = 0;
fdout << "MPFACONNS" << std::endl;
for (Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
{
Cell * r0 = f->BackCell(), *r1 = f->FrontCell();
if(needBuild(&*f))
{
if( f->GetMarker(bnd_markers) && (!bnd_conds.isValid() || !f->HaveData(bnd_conds))) continue;
fdout << f->ReferenceArrayDF(elems)[0]->IntegerDF(unknown_id) << " ";
fdout << f->ReferenceArrayDF(elems)[1]->IntegerDF(unknown_id) << " ";
fdout << f->RealArrayDF(trans)[0]*trans_scale << std::endl;
f->IntegerDF(row) = cnt++;
}
}
fdout << "/" << std::endl;
if( cntsupp )
{
fdout << "SUPPCONNS" << std::endl;
fdout << cntsupp << std::endl;
for (Mesh::iteratorElement f = m->BeginElement(FACE); f != m->EndElement(); ++f)
{
if( f->GetMarker(bnd_markers) && f->HaveData(bnd_conds) )
{
fdout << "--bnd " << ElementTypeName(f->GetElementType()) << "# " << f->LocalID() << std::endl;
fdout << 2 << " "; //boundary type
fdout << f->IntegerDF(row) << " "; //corresponding entry in connection list
// \alpha C + \beta dC/dn = \gamma
if( bnd_conds.isValid() && f->HaveData(bnd_conds) )
{
Storage::real_array bnd = f->RealArray(bnd_conds);
fdout << bnd[0] << " "; // \alpha
fdout << bnd[1] << " "; // \beta
fdout << bnd[2] << " "; // \gamma
}
else
{
fdout << 0.0 << " "; // \alpha
fdout << 1.0 << " "; // \beta
fdout << 0.0 << " "; // \gamma
}
fdout << std::endl;
}
}
fdout << "/" << std::endl;
}
}
void mpfa_o::Update()
{
for (INMOST_DATA_INTEGER_TYPE id = 0; id < m->MaxLocalID(FACE); ++id)
{
Face * f = m->FaceByLocalID(id);
if (f != NULL)
{
if (needBuild(f))
{
Cell * r0 = f->BackCell();
Cell * r1 = f->FrontCell();
if (r0->New() || r1->New()) Build(f);
}
}
}
}
mpfa_o::~mpfa_o()
{
m->DeleteTag(elems);
m->DeleteTag(trans);
}
\ No newline at end of file
#include "discr.h"
#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; }
nmpfa_a::nmpfa_a(Automatizator * aut, Mesh * m, std::string tensor)
:basic_a(aut, m, tensor)
{
}
expr nmpfa_a::Grad(const expr & param) const
{
expr T = tagval(add,0);
expr dp = stencil(trpl, param);
expr dn = stencil(trpr, param);
expr dp2 = ad_abs(dp) + 1.0e-7;
expr dn2 = ad_abs(dn) + 1.0e-7;
expr mp = ad_val(dn2 / (dp2+dn2),0.85);
expr mn = ad_val(dp2 / (dp2+dn2),0.85);
/*
expr dp2 = dp*dp / ad_sqrt(dp*dp + 1.0e-5) + 1.0e-7;
expr dn2 = dn*dn / ad_sqrt(dn*dn + 1.0e-5) + 1.0e-7;
expr mp = ad_val(dn2 / (dn2 + dp2),0.5);
expr mn = ad_val(dp2 / (dn2 + dp2),0.5);
*/
return T*(stencil(upw,param) - stencil(dwn, param)) + dp*mp+dn*mn;
}
expr nmpfa_a::Interp(ElementType etype,const expr & param) const
{
return stencil(intrp,param);
}
void nmpfa_a::Export(std::ostream & fcout, Storage::real trans_scale, Storage::real vol_scale) const
{
discr_basic::Export(fcout,trans_scale,vol_scale);
Tag unknown_id = m->CreateTag("ID",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
Tag row = m->CreateTag("ROW",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
int cnt = 0, cntsupp = 0;
INMOST_DATA_ENUM_TYPE idnum = 0;
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
it->IntegerDF(unknown_id) = idnum++;
for(Mesh::iteratorElement it = m->BeginElement(FACE | EDGE | NODE); it != m->EndElement(); ++it)
if(it->GetMarker(add_markers)) it->IntegerDF(unknown_id) = idnum+cntsupp++;
fcout << std::setprecision(16);
fcout << "NMPFACONNS" << std::endl;
for (Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
{
if (needBuild(&*f))
{
if( f->GetMarker(bnd_markers) && !f->GetMarker(add_markers) ) continue;
cnt++;
}
}
fcout << cnt << std::endl;
cnt = 0;
MIDType marker = m->CreateMarker(), mrkunique = m->CreateMarker();
for (Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
if( f->GetMarker(marker) ) std::cout << "FACE " << f->LocalID() << std::endl;
dynarray<Element *,128> stack;
Tag num = m->CreateTag("TEMP_POS",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
Tag coe = m->CreateTag("TEMP_COE",DATA_REAL,CELL|FACE|EDGE|NODE,NONE,4);
for (Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
{
if (needBuild(&*f))
{
Element * r0 = f->BackCell(), *r1 = f->FrontCell();
int abnP[2] = {0,0};
if( f->GetMarker(bnd_markers) && !f->GetMarker(add_markers) ) continue;
//fcout << "-- face# " << f->LocalID() << " conn# " << cnt;
//if( f->GetMarker(bnd_markers) ) fcout << " is boundary ";
//fcout << std::endl;
{
stack.clear();
stack.push_back(r0);
r0->SetMarker(marker);
r0->IntegerDF(num) = 0;
if( r1 == NULL ) r1 = &*f;
stack.push_back(r1);
r1->SetMarker(marker);
r1->IntegerDF(num) = 1;
Storage::reference_array lrtrp[2];
lrtrp[0] = f->ReferenceArrayDF(trpl_elems);
lrtrp[1] = f->ReferenceArrayDF(trpr_elems);
Storage::real_array lrcoef[2];
lrcoef[0] = f->RealArrayDF(trpl_coefs);
lrcoef[1] = f->RealArrayDF(trpr_coefs);
Storage::real T = f->RealDF(trans);