Skip to content

Multiple sinker model: add option to load inclusion origins from file #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
psanan wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Multiple sinker model: add option to load inclusion origins from file #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7c9b05d
Multiple sinker model: add option to load inclusion origins from file
psanan Apr 12, 2022
43aeb95
Models: add options to control matrix and inclusion densities
psanan Apr 21, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion exSaddle.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ static char help[] = "Solves "
#endif
"-model 6 : single sinker \n"
#if NSD == 3
"-model 7 : lots of sinkers, controlled by -sinker_n \n"
"-model 7 : lots of sinkers, controlled by -sinker_n and -sinker_use_file\n"
#endif
#if defined(LAME)
"-model 8 : one inclusion, fixed base, free elsewhere \n"
Expand Down
92 changes: 70 additions & 22 deletions models.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -655,7 +655,7 @@ Lame' (Elasticity)
PetscErrorCode LameOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu,PetscReal *lambda, PetscReal Fu[],PetscReal Fp[])
{
PetscErrorCode ierr;
static PetscReal opts_mu0,opts_mu1,opts_lambda0,opts_lambda1,opts_rad;
static PetscReal opts_mu0,opts_mu1,opts_lambda0,opts_lambda1,opts_rad,opts_rho0,opts_rho1;
static PetscBool been_here = PETSC_FALSE;
PetscReal mu_qp,lambda_qp,rho_qp;
PetscBool inside = PETSC_FALSE;
Expand All @@ -668,6 +668,8 @@ PetscErrorCode LameOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu
opts_lambda0 = 1.0;
opts_lambda1 = 2.0;
opts_rad = 0.25;
opts_rho0 = 1.0;
opts_rho1 = 2.0;
ierr = PetscOptionsGetReal(NULL,NULL,"-mu0",&opts_mu0,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-mu1",&opts_mu1,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-lambda0",&opts_lambda0,0);CHKERRQ(ierr);
Expand All @@ -683,7 +685,7 @@ PetscErrorCode LameOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu

mu_qp = opts_mu0;
lambda_qp = opts_lambda0;
rho_qp = 1.0;
rho_qp = opts_rho0;

{
#if NSD == 2
Expand All @@ -695,7 +697,7 @@ PetscErrorCode LameOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu
}

if (inside) {
rho_qp = 2.0; /* 2x density inside the inclusion */
rho_qp = opts_rho1;
mu_qp = opts_mu1;
lambda_qp = opts_lambda1;
}
Expand All @@ -722,7 +724,7 @@ PetscErrorCode LameOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu
/******************************************************************************/
PetscErrorCode LameXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu,PetscReal *lambda,PetscReal Fu[],PetscReal Fp[])
{
static PetscReal opts_mu0,opts_mu1,opts_lambda0,opts_lambda1,opts_rad;
static PetscReal opts_mu0,opts_mu1,opts_lambda0,opts_lambda1,opts_rad,opts_rho0,opts_rho1;
static PetscBool been_here = PETSC_FALSE;
static PetscInt opts_numinc;
const PetscReal posx[8] = { 0.27 , 0.6 , 0.7 , 0.2 , 0.85 , 0.4 , 0.16 , 0.55 };
Expand All @@ -744,6 +746,8 @@ PetscErrorCode LameXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu,P
opts_lambda1 = 1.0;
opts_rad = 0.05;
opts_numinc = 3;
opts_rho0 = 1.0;
opts_rho1 = 2.0;

ierr = PetscOptionsGetReal(NULL,NULL,"-mu0",&opts_mu0,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-mu1",&opts_mu1,0);CHKERRQ(ierr);
Expand All @@ -769,7 +773,7 @@ PetscErrorCode LameXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu,P

mu_qp = opts_mu0;
lambda_qp = opts_lambda0;
rho_qp = 1.0;
rho_qp = opts_rho0;

for(i=0; i<opts_numinc; ++i){
#if NSD == 2
Expand All @@ -786,7 +790,7 @@ PetscErrorCode LameXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu,P
if (inside) {
mu_qp = opts_mu1;
lambda_qp = opts_lambda1;
rho_qp = 1.1;
rho_qp = opts_rho1;
}

if (mu) {
Expand All @@ -812,7 +816,7 @@ PetscErrorCode LameXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *mu,P
PetscErrorCode LameHomogeneous_EvaluateCoefficients(PetscReal coor[],PetscReal *mu,PetscReal *lambda, PetscReal Fu[],PetscReal Fp[])
{
PetscErrorCode ierr;
static PetscReal opts_mu0,opts_lambda0;
static PetscReal opts_mu0,opts_lambda0,opts_rho0;
static PetscBool been_here = PETSC_FALSE;
PetscReal mu_qp,lambda_qp,rho_qp;

Expand All @@ -821,6 +825,7 @@ PetscErrorCode LameHomogeneous_EvaluateCoefficients(PetscReal coor[],PetscReal *
PetscPrintf(PETSC_COMM_WORLD,"ModelType: LameHomogeneous\n");
opts_mu0 = 1.0;
opts_lambda0 = 1.0;
opts_rho0 = 1.0;
ierr = PetscOptionsGetReal(NULL,NULL,"-mu0",&opts_mu0,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-lambda0",&opts_lambda0,0);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," params: mu0 %1.4e\n",opts_mu0);CHKERRQ(ierr);
Expand All @@ -830,7 +835,7 @@ PetscErrorCode LameHomogeneous_EvaluateCoefficients(PetscReal coor[],PetscReal *

mu_qp = opts_mu0;
lambda_qp = opts_lambda0;
rho_qp = 1.0;
rho_qp = opts_rho0;

if (lambda) {
*lambda = lambda_qp;
Expand Down Expand Up @@ -935,7 +940,7 @@ PetscErrorCode StokesSolCx_EvaluateCoefficients(PetscReal coor[],PetscReal *eta,
/*****************************************************************************/
PetscErrorCode StokesThreeSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *eta,PetscReal Fu[],PetscReal Fp[])
{
static PetscReal opts_eta0,opts_eta1,opts_rad;
static PetscReal opts_eta0,opts_eta1,opts_rad,opts_rho0,opts_rho1;
static PetscBool been_here = PETSC_FALSE;
PetscReal eta_qp,rho_qp;
PetscBool inside = PETSC_FALSE;
Expand All @@ -947,6 +952,8 @@ PetscErrorCode StokesThreeSinker_EvaluateCoefficients(PetscReal coor[],PetscReal
opts_eta0 = 1.0;
opts_eta1 = 1.0;
opts_rad = 0.1;
opts_rho0 = 1.0;
opts_rho1 = 1.1;

ierr = PetscOptionsGetReal(NULL,NULL,"-eta0",&opts_eta0,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-eta1",&opts_eta1,0);CHKERRQ(ierr);
Expand All @@ -958,7 +965,7 @@ PetscErrorCode StokesThreeSinker_EvaluateCoefficients(PetscReal coor[],PetscReal
}

eta_qp = opts_eta0;
rho_qp = 1.0;
rho_qp = opts_rho0;

{
PetscReal sep2;
Expand All @@ -985,7 +992,7 @@ PetscErrorCode StokesThreeSinker_EvaluateCoefficients(PetscReal coor[],PetscReal

if (inside) {
eta_qp = opts_eta1;
rho_qp = 1.1;
rho_qp = opts_rho1;
}

if (eta) {
Expand All @@ -1006,7 +1013,7 @@ PetscErrorCode StokesThreeSinker_EvaluateCoefficients(PetscReal coor[],PetscReal
/*****************************************************************************/
PetscErrorCode StokesXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *eta,PetscReal Fu[],PetscReal Fp[])
{
static PetscReal opts_eta0,opts_eta1,opts_rad;
static PetscReal opts_eta0,opts_eta1,opts_rad,opts_rho0,opts_rho1;
static PetscBool been_here = PETSC_FALSE;
static PetscInt opts_numinc;
const PetscReal posx[8] = { 0.27 , 0.6 , 0.7 , 0.2 , 0.85 , 0.4 , 0.16 , 0.55 };
Expand All @@ -1026,6 +1033,8 @@ PetscErrorCode StokesXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *et
opts_eta1 = 1.0;
opts_rad = 0.05;
opts_numinc = 3;
opts_rho0 = 1.0;
opts_rho1 = 1.1;

ierr = PetscOptionsGetReal(NULL,NULL,"-eta0",&opts_eta0,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-eta1",&opts_eta1,0);CHKERRQ(ierr);
Expand All @@ -1046,7 +1055,7 @@ PetscErrorCode StokesXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *et
}

eta_qp = opts_eta0;
rho_qp = 1.0;
rho_qp = opts_rho0;

for(i=0; i<opts_numinc; ++i){
#if NSD == 2
Expand All @@ -1062,7 +1071,7 @@ PetscErrorCode StokesXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *et

if (inside) {
eta_qp = opts_eta1;
rho_qp = 1.1;
rho_qp = opts_rho1;
}

if (eta) {
Expand All @@ -1084,7 +1093,7 @@ PetscErrorCode StokesXSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *et
/*****************************************************************************/
PetscErrorCode StokesOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *eta,PetscReal Fu[],PetscReal Fp[])
{
static PetscReal opts_eta0,opts_eta1,opts_rad,opts_x,opts_y;
static PetscReal opts_eta0,opts_eta1,opts_rad,opts_x,opts_y,opts_rho0,opts_rho1;
#if NSD==3
static PetscReal opts_z;
#endif
Expand All @@ -1104,6 +1113,8 @@ PetscErrorCode StokesOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *
#if NSD==3
opts_z = 0.5;
#endif
opts_rho0 = 1.0;
opts_rho1 = 1.1;

ierr = PetscOptionsGetReal(NULL,NULL,"-eta0",&opts_eta0,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-eta1",&opts_eta1,0);CHKERRQ(ierr);
Expand All @@ -1125,7 +1136,7 @@ PetscErrorCode StokesOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *
}

eta_qp = opts_eta0;
rho_qp = 1.0;
rho_qp = opts_rho0;

{
PetscReal sep2;
Expand All @@ -1139,7 +1150,7 @@ PetscErrorCode StokesOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *

if (inside) {
eta_qp = opts_eta1;
rho_qp = 1.1;
rho_qp = opts_rho1;
}

if (eta) {
Expand All @@ -1160,6 +1171,34 @@ PetscErrorCode StokesOneSinker_EvaluateCoefficients(PetscReal coor[],PetscReal *

/*****************************************************************************/
#if NSD == 3
/*
(3d only)
Loads inclusion origins from a simple text file, with three floating point numbers per line.
*/
PetscErrorCode LoadInclusionOrigins(const char* filename,PetscInt ninclusions,PetscReal **_pos)
{
PetscReal *pos;
FILE *fp;

PetscPrintf(PETSC_COMM_WORLD,"# LoadInclusionOrigins\n");
PetscPrintf(PETSC_COMM_WORLD,"Loading %D inclusion origin(s) from file: %s\n",ninclusions,filename);

PetscMalloc1(NSD*ninclusions,&pos);
if ((fp = fopen(filename,"r")) == NULL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Could not open inclusions file %s",filename);
for (PetscInt i=0; i<ninclusions; ++i) {
double xp, yp, zp;

if (fscanf(fp,"%lf %lf %lf",&xp,&yp,&zp) != 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Could not read coordinates for inclusion %D. Does the file have enough lines?",i+1);
PetscPrintf(PETSC_COMM_WORLD,"%lf %lf %lf\n",xp,yp,zp);
/* Here, one could check that the points make sense with the specified radius, etc. */
pos[NSD*i+0] = xp;
pos[NSD*i+1] = yp;
pos[NSD*i+2] = zp;
}
*_pos = pos;
PetscFunctionReturn(0);
}

/*
(3d only)
Generates location of spherical inclusions in the domain [0,Lx]x[0,Ly]x[0,Lz] assuming the
Expand Down Expand Up @@ -1267,9 +1306,10 @@ PetscErrorCode GenerateInclusionOrigins(PetscInt ninclusions,PetscReal rmax,Pets
PetscErrorCode SinkerPtatin_EvaluateCoefficients(PetscReal coor[],PetscReal *eta,PetscReal Fu[],PetscReal Fp[])
{
PetscErrorCode ierr;
static PetscReal opts_eta0,opts_eta1,opts_rad,*centroidpos;
static PetscBool been_here = PETSC_FALSE;
static PetscReal opts_eta0,opts_eta1,opts_rad,*centroidpos,opts_rho0,opts_rho1;
static PetscBool opts_use_file = PETSC_FALSE, been_here = PETSC_FALSE;
static PetscInt opts_numinc;
static char opts_inclusions_filename[PETSC_MAX_PATH_LEN];
PetscReal eta_qp,rho_qp;
PetscBool inside;
PetscInt k;
Expand All @@ -1281,21 +1321,29 @@ PetscErrorCode SinkerPtatin_EvaluateCoefficients(PetscReal coor[],PetscReal *eta
opts_eta1 = 1.1;
opts_rad = 0.05;
opts_numinc = 3;
opts_rho0 = 1.0;
opts_rho1 = 1.1;

ierr = PetscOptionsGetReal(NULL,NULL,"-eta0",&opts_eta0,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-eta1",&opts_eta1,0);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_rad,0);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-sinker_n",&opts_numinc,0);CHKERRQ(ierr);
ierr = PetscOptionsGetString(NULL,NULL,"-sinker_use_file",opts_inclusions_filename,PETSC_MAX_PATH_LEN-1,&opts_use_file);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD," params: eta0 %1.4e\n",opts_eta0);
PetscPrintf(PETSC_COMM_WORLD," params: eta1 %1.4e\n",opts_eta1);

ierr = GenerateInclusionOrigins(opts_numinc,opts_rad,1.0,1.0,1.0,1.5,1.5,&centroidpos);CHKERRQ(ierr);
if (opts_use_file) {
ierr = LoadInclusionOrigins(opts_inclusions_filename,opts_numinc,&centroidpos);CHKERRQ(ierr);
} else {
ierr = GenerateInclusionOrigins(opts_numinc,opts_rad,1.0,1.0,1.0,1.5,1.5,&centroidpos);CHKERRQ(ierr);
/* Here, one could consider dumping out the inclusion origins to a text file to use later */
}

been_here = PETSC_TRUE;
}

eta_qp = opts_eta0;
rho_qp = 1.0;
rho_qp = opts_rho0;

inside = PETSC_FALSE;
for (k=0; k<opts_numinc; ++k) {
Expand All @@ -1316,7 +1364,7 @@ PetscErrorCode SinkerPtatin_EvaluateCoefficients(PetscReal coor[],PetscReal *eta

if (inside) {
eta_qp = opts_eta1;
rho_qp = 1.1;
rho_qp = opts_rho1;
}

if (eta) {
Expand Down