OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-22 15:33:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Split mim_ip_simple() up into two co-operating procedures,

Browse Source 
- mim_ip_span_nullspace()
        Constructs orthogonal basis for mimetic IP null space.
    - mim_ip_linpress_exact()
        Adds term to impose exactness for linear pressure fields.

  Re-implement mim_ip_simple() in terms of these procedures.
This commit is contained in:
Bård Skaflestad 2010-06-29 21:18:23 +00:00
parent 9e1bbb200b
commit da327ba410
2 changed files with 50 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

95
mimetic.c
View File

@ -1,71 +1,41 @@
#include <stddef.h>
#ifndef MAT_SIZE_T
#define MAT_SIZE_T size_t
#define MAT_SIZE_T int
#endif
#include "blas_lapack.h"
#include "mimetic.h"
/* ------------------------------------------------------------------ */
static double
trace(size_t n, double *A);
/* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
static void
mim_ip_simple_fill_null(size_t nf, size_t d, double *C, double *A,
double *Binv, double *work, MAT_SIZE_T lwork);
/* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
void
mim_ip_simple(MAT_SIZE_T nf, MAT_SIZE_T d, double vol,
mim_ip_simple(int nf, int d, double vol,
double *A, double *K, double *C, double *N,
double *Binv, double* work, MAT_SIZE_T lwork)
double *Binv, double *work, int lwork)
/* ------------------------------------------------------------------ */
{
MAT_SIZE_T m, n, k, ld1, ld2, ldBinv;
double a1, a2;
mim_ip_simple_fill_null(nf, d, C, A, Binv, work, lwork);
/* Step 4) C <- N*K */
m = nf ; n = nf ; k = d;
ld1 = nf ; ld2 = d ;
a1 = 1.0; a2 = 0.0;
dgemm_("No Transpose", "No Transpose", &m, &n, &k,
&a1, N, &ld1, K, &ld1, &a2, C, &ld1);
/* Step 5) Binv <- (N*K*N' + t*A*(I-Q*Q')*A) / vol */
a1 = 1.0 / vol ;
a2 = 6.0 * trace(nf, K) / (d * vol);
ldBinv = nf;
dgemm_("No Transpose", "Transpose", &m, &m, &n,
&a1, C, &ld1, N, &ld1, &a2, Binv, &ldBinv);
mim_ip_span_nullspace(nf, d, C, A, Binv, work, lwork);
mim_ip_linpress_exact(nf, d, vol, N, K, Binv, C);
}
/* ------------------------------------------------------------------ */
static void
mim_ip_simple_fill_null(size_t nf, size_t d, double *C, double *A,
double *Binv,
double *work, MAT_SIZE_T lwork)
void
mim_ip_span_nullspace(int nf, int d, double *C, double *A,
double *X, double *work, int nwork)
/* ------------------------------------------------------------------ */
{
MAT_SIZE_T m, n, k, ld, info;
MAT_SIZE_T m, n, k, ld, info, lwork;
size_t i, j;
int i, j;
double a1, a2;
double tau[3] = { 0.0 }; /* No more than 3 spatial dimensions */
/* Step 1) Binv <- I */
for (i = 0; i < nf*nf; i++) { Binv[i ] = 0.0; }
for (i = 0; i < nf ; i++) { Binv[i * (nf + 1)] = 1.0; }
/* Step 1) X <- I */
for (i = 0; i < nf*nf; i++) { X[i ] = 0.0; }
for (i = 0; i < nf ; i++) { X[i * (nf + 1)] = 1.0; }
/* Step 2) C <- orth(A * C) */
for (j = 0; j < d; j++) {
@ -74,40 +44,55 @@ mim_ip_simple_fill_null(size_t nf, size_t d, double *C, double *A,
}
}
m = nf; n = d; ld = nf; k = d;
m = nf; n = d; ld = nf; k = d; lwork = nwork;
dgeqrf_(&m, &n, C, &ld, tau, work, &lwork, &info);
dorgqr_(&m, &n, &k, C, &ld, tau, work, &lwork, &info);
/* Step 3) Binv <- A * (Binv - C*C') * A */
/* Step 3) X <- A * (X - C*C') * A */
a1 = -1.0; a2 = 1.0;
dsyrk_("Upper Triangular", "No Transpose",
&m, &n, &a1, C, &ld, &a2, Binv, &ld);
&m, &n, &a1, C, &ld, &a2, X, &ld);
for (j = 0; j < nf; j++) {
for (i = 0; i <= j; i++) {
Binv[i + j*nf] *= A[i] * A[j];
X[i + j*nf] *= A[i] * A[j];
}
}
/* Account for DSYRK only assigning upper triangular part. */
for (j = 0; j < nf; j++) {
for (i = j + 1; i < nf; i++) {
Binv[i + j*nf] = Binv[j + i*nf];
X[i + j*nf] = X[j + i*nf];
}
}
}
/* ------------------------------------------------------------------ */
static double
trace(size_t n, double *A)
void
mim_ip_linpress_exact(int nf, int d, double vol, double *N, double *K,
double *Binv, double *T)
/* ------------------------------------------------------------------ */
{
size_t i;
double t = 0.0;
MAT_SIZE_T m, n, k, ld1, ld2, ldBinv;
int i;
double a1, a2, t;
for (i = 0; i < n; i++) {
t += A[i + i*n];
t = 0.0;
for (i = 0; i < nf; i++) {
t += K[i + i*nf];
}
return t;
/* Step 4) T <- N*K */
m = nf ; n = nf ; k = d;
ld1 = nf ; ld2 = d ;
a1 = 1.0; a2 = 0.0;
dgemm_("No Transpose", "No Transpose", &m, &n, &k,
&a1, N, &ld1, K, &ld1, &a2, T, &ld1);
/* Step 5) Binv <- (N*K*N' + t*X) / vol */
a1 = 1.0 / vol ;
a2 = 6.0 * t / (d * vol);
ldBinv = nf;
dgemm_("No Transpose", "Transpose", &m, &m, &n,
&a1, T, &ld1, N, &ld1, &a2, Binv, &ldBinv);
}

14
mimetic.h
View File

@ -1,10 +1,16 @@
#ifndef MIMETIC_H_INCLUDED
#define MIMETIC_H_INCLUDED
void mim_ip_simple(MAT_SIZE_T nf, MAT_SIZE_T d, double v,
double *A, double *K, double *C, double *N,
double *Binv,
double *work, MAT_SIZE_T lwork);
void mim_ip_span_nullspace(int nf, int d, double *C, double *A,
double *X, double *work, int lwork);
void mim_ip_linpress_exact(int nf, int d, double vol, double *N,
double *K, double *Binv, double *T);
void mim_ip_simple(int nf, int d, double v, double *A, double *K,
double *C, double *N, double *Binv,
double *work, int lwork);
#endif /* MIMETIC_H_INCLUDED */