OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Updated files from MRST repository.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2012-01-20 13:09:13 +01:00
parent 7d6fb03142
commit db8681b9ff
4 changed files with 77 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files

80
opm/core/transport/reorder/reordersequence.c
View File

@@ -2,34 +2,42 @@
#include <stdio.h>
#include <stdlib.h>
#ifdef MATLAB_MEX_FILE
#include "grid.h"
#include "reordersequence.h"
#include "tarjan.h"
#else
#include <opm/core/grid.h>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/transport/reorder/tarjan.h>
#endif
/* Construct adjacency matrix of upwind graph wrt flux. Column
indices are not sorted. */
static void
make_upwind_graph(int nc, int *cellfaces, int *faceptr, int *face2cell,
static void
make_upwind_graph(int nc, int *cellfaces, int *faceptr, int *face2cell,
const double *flux, int *ia, int *ja, int *work)
{
/* Using topology (conn, cptr), and direction, construct adjacency
matrix of graph. */
int i, j, p, f, positive_sign;
int i, j, p, f, positive_sign, boundaryface;
double theflux;
/* For each face, store upwind cell in work array */
for (i=0; i<nc; ++i) {
for (j=faceptr[i]; j<faceptr[i+1]; ++j) {
for (i=0; i<nc; ++i)
{
for (j=faceptr[i]; j<faceptr[i+1]; ++j)
{
f = cellfaces[j];
positive_sign = (i == face2cell[2*f]);
theflux = positive_sign ? flux[f] : -flux[f];
if ( theflux > 0 ) {
/* flux from cell i over face f */
work[f] = i;
if ( theflux > 0 )
{
/* i is upwind cell for face f */
work[f] = i;
}
}
}
@@ -37,14 +45,25 @@ make_upwind_graph(int nc, int *cellfaces, int *faceptr, int *face2cell,
/* Fill ia and ja */
p = 0;
ia[0] = p;
for (i=0; i<nc; ++i) {
for (j=faceptr[i]; j<faceptr[i+1]; ++j) {
for (i=0; i<nc; ++i)
{
for (j=faceptr[i]; j<faceptr[i+1]; ++j)
{
f = cellfaces[j];
boundaryface = (face2cell[2*f+0] == -1) ||
(face2cell[2*f+1] == -1);
if ( boundaryface )
{
continue;
}
positive_sign = (i == face2cell[2*f]);
theflux = positive_sign ? flux[f] : -flux[f];
if ( theflux < 0) {
if ( theflux < 0)
{
ja[p++] = work[f];
}
}
@@ -53,45 +72,46 @@ make_upwind_graph(int nc, int *cellfaces, int *faceptr, int *face2cell,
}
}
static void
compute_reorder_sequence(int nc, int nf, int *cellfaces, int *faceptr, int *face2cell,
static void
compute_reorder_sequence(int nc, int nf, int *cellfaces, int *faceptr, int *face2cell,
const double *flux, int *sequence, int *components, int *ncomponents)
{
int *ia, *ja;
int *work;
int sz = nf;
if (nf < 3*nc) {
if (nf < 3*nc)
{
sz = 3*nc;
}
work = malloc( sz * sizeof *work);
ia = malloc((nc+1) * sizeof *ia);
ja = malloc( nf * sizeof *ja); /* A bit too much... */
make_upwind_graph(nc, cellfaces, faceptr, face2cell,
make_upwind_graph(nc, cellfaces, faceptr, face2cell,
flux, ia, ja, work);
tarjan (nc, ia, ja, sequence, components, ncomponents, work);
free(ja);
free(ia);
free(work);
}
void compute_sequence(struct UnstructuredGrid *grid, const double *flux,
int *sequence,
void compute_sequence(struct UnstructuredGrid *grid, const double *flux,
int *sequence,
int *components, int *ncomponents)
{
compute_reorder_sequence(grid->number_of_cells,
grid->number_of_faces,
grid->cell_faces,
grid->cell_facepos,
compute_reorder_sequence(grid->number_of_cells,
grid->number_of_faces,
grid->cell_faces,
grid->cell_facepos,
grid->face_cells,
flux,
sequence,
components,
flux,
sequence,
components,
ncomponents);
}
@@ -99,5 +119,3 @@ void compute_sequence(struct UnstructuredGrid *grid, const double *flux,
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

2
opm/core/transport/reorder/reordersequence.h
View File

@@ -3,7 +3,7 @@
#define REORDERSEQUENCE_H_INCLUDED
struct UnstructuredGrid;
void compute_sequence(struct UnstructuredGrid *grid, const double *flux,
void compute_sequence(struct UnstructuredGrid *grid, const double *flux,
int *sequence, int *components, int *ncomponents);
#endif

84
opm/core/transport/reorder/tarjan.c
View File

@@ -2,13 +2,17 @@
#include <string.h>
#include <assert.h>
#ifdef MATLAB_MEX_FILE
#include "tarjan.h"
#else
#include <opm/core/transport/reorder/tarjan.h>
#endif
static int min(int a, int b){ return a<b? a : b;}
static int min(int a, int b){ return a < b? a : b;}
/* Improved (or obfuscated) version uses less memory
*
@@ -36,7 +40,7 @@ tarjan (int size, int *ia, int *ja, int *P, int *Q, int *ncomp,
int *link = (int *) time + size;
int *status = (int*) link + size; /* dual usage... */
(void) cbottom;
memset(work, 0, 3*size * sizeof *work);
memset(P, 0, size * sizeof *P );
@@ -50,27 +54,27 @@ tarjan (int size, int *ia, int *ja, int *P, int *Q, int *ncomp,
*ncomp = 0;
*Q++ = pos;
seed = 0;
while (seed < size)
{
{
if (status[seed] == DONE)
{
++seed;
continue;
}
*stack-- = seed; /* push seed */
t = 0;
while ( stack != bottom )
{
c = *(stack+1); /* peek c */
assert(status[c] != DONE);
assert(status[c] >= -2);
if (status[c] == REMAINING)
{
status[c] = ia[c+1]-ia[c]; /* number of descendants of c */
@@ -81,11 +85,11 @@ tarjan (int size, int *ia, int *ja, int *P, int *Q, int *ncomp,
/* if all descendants are processed */
if (status[c] == 0)
if (status[c] == 0)
{
/* if c is root of strong component */
if (link[c] == time[c])
if (link[c] == time[c])
{
do
{
@@ -94,16 +98,16 @@ tarjan (int size, int *ia, int *ja, int *P, int *Q, int *ncomp,
v = *++cstack; /* pop strong component stack */
status[v] = DONE;
P[pos++] = v; /* store vertex in P */
}
}
while ( v != c );
*Q++ = pos; /* store end point of component */
*Q++ = pos; /* store end point of component */
++*ncomp;
}
++stack; /* pop c */
if (stack != bottom)
if (stack != bottom)
{
link[*(stack+1)] = min(link[*(stack+1)], link[c]);
}
@@ -112,24 +116,24 @@ tarjan (int size, int *ia, int *ja, int *P, int *Q, int *ncomp,
/* if there are more descendants to consider */
else
else
{
assert(status[c] > 0);
child = ja[ia[c] + status[c]-1];
--status[c]; /* Decrement descendant count of c*/
if (status[child] == REMAINING)
if (status[child] == REMAINING)
{
*stack-- = child; /* push child */
}
else if (status[child] >= 0)
}
else if (status[child] >= 0)
{
link[c] = min(link[c], time[child]);
}
else
}
else
{
assert(status[child] = DONE);
}
@@ -139,42 +143,6 @@ tarjan (int size, int *ia, int *ja, int *P, int *Q, int *ncomp,
}
}
#ifdef TEST_TARJAN
#include <stdio.h>
#include <stdlib.h>
/* gcc -DTEST_TARJAN tarjan.c */
int main()
{
int size = 5;
int ia[] = {0, 1, 2, 4, 5, 5};
int ja[] = {1, 2, 0, 3, 4};
int *work = malloc(3*size*sizeof *work);
int number_of_components;
int *permutation = malloc(size * sizeof *permutation);
int *ptr = malloc(size * sizeof *ptr);
int i,j;
tarjan (size, ia, ja, permutation, ptr, &number_of_components, work);
fprintf(stderr, "Number of components: %d\n", number_of_components);
for (i=0; i<number_of_components; ++i)
{
for (j=ptr[i]; j<ptr[i+1]; ++j)
{
fprintf(stderr, "%d %d\n", i, permutation[j]);
}
}
free(work);
free(permutation);
free(ptr);
return 0;
}
#endif
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

2
opm/core/transport/reorder/tarjan.h
View File

@@ -6,7 +6,7 @@
extern "C" {
#endif
void tarjan (int size, int *ia, int *ja, int *rowP, int *P,
void tarjan (int size, int *ia, int *ja, int *rowP, int *P,
int *ncomp, int *work);
#ifdef __cplusplus