Skip to content

Commit

Permalink
nn
Browse files Browse the repository at this point in the history
  • Loading branch information
@fegomezl
fegomezl committed May 30, 2022
1 parent cb3d490 commit c28e119
Show file tree
Hide file tree
Showing 9 changed files with 76 additions and 130 deletions.
43 changes: 18 additions & 25 deletions Brinicle/code/assemble.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,13 @@

//Usefull position functions
double r_f(const Vector &x);
double r_inv_f(const Vector &x);
void zero_f(const Vector &x, Vector &f);
void r_inv_hat_f(const Vector &x, Vector &f);
void rot_f(const Vector &x, DenseMatrix &f);

//Physical properties (in T,S)
double FusionPoint(const double T, const double S);
double FusionPoint(const double S);
double RelativeTemperature(const double T, const double S);
double Phase(const double T, const double S);
double HeatDiffusivity(const double T, const double S);
double SaltDiffusivity(const double T, const double S);
Expand All @@ -21,34 +21,30 @@ void Artic_sea::assemble_system(){
//Define fields
temperature = new ParGridFunction(fespace_H1);
salinity = new ParGridFunction(fespace_H1);
relative_temperature = new ParGridFunction(fespace_H1);
phase = new ParGridFunction(fespace_H1);
vorticity = new ParGridFunction(fespace_H1);
stream = new ParGridFunction(fespace_H1);
velocity = new ParGridFunction(fespace_ND);
rvelocity = new ParGridFunction(fespace_ND);

rVelocity = new HypreParVector(fespace_ND);
Velocity = new HypreParVector(fespace_ND);

//Initialize operators
transport_oper = new Transport_Operator(config, *fespace_H1, *fespace_ND, dim, pmesh->bdr_attributes.Max(), block_offsets_H1, X);
flow_oper = new Flow_Operator(config, *fespace_H1, *fespace_ND, dim, pmesh->bdr_attributes.Max(), block_offsets_H1);
transport_oper = new Transport_Operator(config, *fespace_H1, dim, pmesh->bdr_attributes.Max(), block_offsets_H1, X);
flow_oper = new Flow_Operator(config, *fespace_H1, dim, pmesh->bdr_attributes.Max(), block_offsets_H1);

//Solve initial velocity field
flow_oper->SetParameters(X);
flow_oper->Solve(Y, *Velocity, *rVelocity);
flow_oper->Solve(Y);

//Set initial state
temperature->Distribute(X.GetBlock(0));
salinity->Distribute(X.GetBlock(1));
vorticity->Distribute(Y.GetBlock(0));
stream->Distribute(Y.GetBlock(1));
velocity->Distribute(Velocity);
rvelocity->Distribute(rVelocity);

//Calculate phases
for (int ii = 0; ii < phase->Size(); ii++)
(*phase)(ii) = FusionPoint((*temperature)(ii), (*salinity)(ii));
for (int ii = 0; ii < phase->Size(); ii++){
(*relative_temperature)(ii) = RelativeTemperature((*temperature)(ii), (*salinity)(ii));
(*phase)(ii) = Phase((*temperature)(ii), (*salinity)(ii));
}

//Set the ODE solver type
arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::IMPLICIT);
Expand All @@ -65,11 +61,10 @@ void Artic_sea::assemble_system(){
paraview_out->SetLevelsOfDetail(config.order);
paraview_out->RegisterField("Temperature", temperature);
paraview_out->RegisterField("Salinity", salinity);
paraview_out->RegisterField("RelativeTemperature", relative_temperature);
paraview_out->RegisterField("Phase", phase);
paraview_out->RegisterField("Vorticity", vorticity);
paraview_out->RegisterField("Stream", stream);
paraview_out->RegisterField("Velocity", velocity);
paraview_out->RegisterField("rVelocity", rvelocity);
paraview_out->SetCycle(vis_print);
paraview_out->SetTime(t*t_ref);
paraview_out->Save();
Expand Down Expand Up @@ -101,11 +96,6 @@ double r_f(const Vector &x){
return x(0);
}

//Function for 1/r
double r_inv_f(const Vector &x){
return pow(x(0), -1);
}

//Function for 0 (vector)
void zero_f(const Vector &x, Vector &f){
f(0) = 0.;
Expand All @@ -126,15 +116,18 @@ void rot_f(const Vector &x, DenseMatrix &f){
}

//Fusion temperature at a given salinity
double FusionPoint(const double T, const double S){
double T_ = T + ZeroTemperature;
double FusionPoint(const double S){
double S_ = S + ZeroSalinity;
return (1-2*signbit(T_ref))*(T_ - (constants.FusionPoint_a*(S_) + constants.FusionPoint_b*pow(S_, 3)));
return constants.FusionPoint_a*(S_) + constants.FusionPoint_b*pow(S_, 3) - ZeroTemperature;
}

double RelativeTemperature(const double T, const double S){
return (1-2*signbit(T_ref))*(T - FusionPoint(S));
}

//Phase indicator (1 for liquid and 0 for solid)
double Phase(const double T, const double S){
return 0.5*(1+tanh(5*EpsilonInv*FusionPoint(T, S)));
return 0.5*(1+tanh(5*EpsilonInv*RelativeTemperature(T, S)));
}

//Coefficient for the diffusion term in the temperature equation
Expand Down
19 changes: 4 additions & 15 deletions Brinicle/code/flow_assemble.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,31 +4,25 @@
double boundary_stream_f(const Vector &x);
double boundary_stream_in_f(const Vector &x);

Flow_Operator::Flow_Operator(Config config, ParFiniteElementSpace &fespace_H1, ParFiniteElementSpace &fespace_ND, int dim, int attributes, Array<int> block_offsets_H1):
Flow_Operator::Flow_Operator(Config config, ParFiniteElementSpace &fespace_H1, int dim, int attributes, Array<int> block_offsets_H1):
config(config),
fespace_H1(fespace_H1),
block_offsets_H1(block_offsets_H1),
ess_bdr_0(attributes), ess_bdr_1(attributes),
ess_bdr_in(attributes),
ess_bdr_closed_down(attributes), ess_bdr_closed_up(attributes),
vorticity(&fespace_H1), stream(&fespace_H1),
stream_gradient(&fespace_ND),
velocity(&fespace_ND), rvelocity(&fespace_ND),
temperature(&fespace_H1), salinity(&fespace_H1),
density(&fespace_H1), density_dr(&fespace_H1),
impermeability(&fespace_H1),
B(block_offsets_H1),
B0(NULL), B1(NULL),
A00(NULL), A01(NULL), A10(NULL), A11(NULL),
coeff_r(r_f), coeff_r_inv(r_inv_f),
coeff_r_inv_hat(dim, r_inv_hat_f),
coeff_rot(dim, rot_f),
coeff_vorticity(0.),
coeff_stream(boundary_stream_f),
coeff_r(r_f), coeff_r_inv_hat(dim, r_inv_hat_f),
coeff_vorticity(0.), coeff_stream(boundary_stream_f),
coeff_stream_in(boundary_stream_in_f),
coeff_stream_closed_down(0.),
coeff_stream_closed_up(1.),
gradient(&fespace_H1, &fespace_ND)
coeff_stream_closed_up(1.)
{
/****
* Define essential boundary conditions
Expand Down Expand Up @@ -108,11 +102,6 @@ Flow_Operator::Flow_Operator(Config config, ParFiniteElementSpace &fespace_H1, P
A01 = a01.ParallelAssemble();

B0 = b0.ParallelAssemble();

//Create gradient interpolator
gradient.AddDomainIntegrator(new GradientInterpolator);
gradient.Assemble();
gradient.Finalize();
}

//Boundary condition for stream
Expand Down
14 changes: 1 addition & 13 deletions Brinicle/code/flow_solve.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#include "header.h"

//Solution of the current system
void Flow_Operator::Solve(BlockVector &Y, Vector &Velocity, Vector &rVelocity){
void Flow_Operator::Solve(BlockVector &Y){

//Create the complete bilinear operator:
//
Expand Down Expand Up @@ -35,17 +35,5 @@ void Flow_Operator::Solve(BlockVector &Y, Vector &Velocity, Vector &rVelocity){
vorticity.Distribute(Y.GetBlock(0));
stream.Distribute(Y.GetBlock(1));

//Calculate velocity field
gradient.Mult(stream, stream_gradient);
VectorGridFunctionCoefficient coeff_stream_gradient(&stream_gradient);
MatrixVectorProductCoefficient coeff_rV(coeff_rot, coeff_stream_gradient);
ScalarVectorProductCoefficient coeff_V(coeff_r_inv, coeff_rV);
velocity.ProjectDiscCoefficient(coeff_V, GridFunction::ARITHMETIC);
rvelocity.ProjectDiscCoefficient(coeff_rV, GridFunction::ARITHMETIC);

//Export flow information
velocity.ParallelAverage(Velocity);
rvelocity.ParallelAverage(rVelocity);

delete H;
}
13 changes: 4 additions & 9 deletions Brinicle/code/grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,10 +37,6 @@ void Artic_sea::make_grid(const char *mesh_file){
fespace_H1 = new ParFiniteElementSpace(pmesh, fec_H1);
size_H1 = fespace_H1->GlobalTrueVSize();

fec_ND = new ND_FECollection(config.order, dim);
fespace_ND = new ParFiniteElementSpace(pmesh, fec_ND);
size_ND = fespace_ND->GlobalTrueVSize();

//Create the block offsets
block_offsets_H1[0] = 0;
block_offsets_H1[1] = fespace_H1->TrueVSize();
Expand All @@ -51,10 +47,9 @@ void Artic_sea::make_grid(const char *mesh_file){
Y.Update(block_offsets_H1); Y = 0.;

if (config.master)
cout << "\nSize (H1): " << size_H1 << "\n"
<< "Size (ND): " << size_ND << "\n"
<< "Mesh Size: " << h_min << " (" << h_min*L_ref << " mm)\n"
<< "Serial refinements: " << serial_refinements << "\n"
cout << "\nSerial refinements: " << serial_refinements << "\n"
<< "Parallel refinements: " << config.refinements - serial_refinements << "\n"
<< "Total refinements: " << config.refinements << "\n\n";
<< "Total refinements: " << config.refinements << "\n"
<< "Size (H1): " << size_H1 << "\n"
<< "Mesh Size: " << h_min << " (" << h_min*L_ref << " mm)\n\n";
}
39 changes: 13 additions & 26 deletions Brinicle/code/header.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,10 +118,10 @@ struct Config{
class Transport_Operator : public TimeDependentOperator{
public:
//Initialization of the solver
Transport_Operator(Config config, ParFiniteElementSpace &fespace_H1, ParFiniteElementSpace &fespace_ND, int dim, int attributes, Array<int> block_offsets_H1, BlockVector &X);
Transport_Operator(Config config, ParFiniteElementSpace &fespace_H1, int dim, int attributes, Array<int> block_offsets_H1, BlockVector &X);

//Update of the solver on each iteration
void SetParameters(const BlockVector &X, const Vector &rVelocity);
void SetParameters(const BlockVector &X, const BlockVector &Y);

//Time-evolving functions
virtual void Mult(const Vector &X, Vector &dX_dt) const;
Expand All @@ -143,20 +143,22 @@ class Transport_Operator : public TimeDependentOperator{
Array<int> ess_bdr_0, ess_bdr_1;

//Auxiliar grid functions
ParGridFunction temperature, salinity, phase;
ParGridFunction rvelocity;
ParGridFunction temperature, salinity;
ParGridFunction relative_temperature, phase;
ParGridFunction stream;
ParGridFunction heat_diffusivity;
ParGridFunction salt_diffusivity;

//Coefficients
FunctionCoefficient coeff_r;
VectorFunctionCoefficient coeff_zero;
MatrixFunctionCoefficient coeff_rot;

ProductCoefficient coeff_rM;
ProductCoefficient coeff_rD0;
ProductCoefficient coeff_rD1;

VectorGridFunctionCoefficient coeff_rV;
MatrixVectorProductCoefficient coeff_rV;
ScalarVectorProductCoefficient coeff_rMV;

InnerProductCoefficient coeff_dPdT;
Expand Down Expand Up @@ -185,13 +187,13 @@ class Transport_Operator : public TimeDependentOperator{
class Flow_Operator{
public:
//Initialization of the solver
Flow_Operator(Config config, ParFiniteElementSpace &fespace_H1, ParFiniteElementSpace &fespace_ND, int dim, int attributes, Array<int> block_offsets_H1);
Flow_Operator(Config config, ParFiniteElementSpace &fespace_H1, int dim, int attributes, Array<int> block_offsets_H1);

//Update of the solver on each iteration
void SetParameters(const BlockVector &X);

//Solution of the current system
void Solve(BlockVector &Y, Vector &Velocity, Vector &rVelocity);
void Solve(BlockVector &Y);

~Flow_Operator();
protected:
Expand All @@ -212,10 +214,6 @@ class Flow_Operator{
//Auxiliar grid functions
ParGridFunction vorticity;
ParGridFunction stream;
ParGridFunction stream_gradient;
ParGridFunction velocity;
ParGridFunction rvelocity;

ParGridFunction temperature;
ParGridFunction salinity;
ParGridFunction density;
Expand All @@ -233,18 +231,13 @@ class Flow_Operator{

//Coefficients
FunctionCoefficient coeff_r;
FunctionCoefficient coeff_r_inv;
VectorFunctionCoefficient coeff_r_inv_hat;
MatrixFunctionCoefficient coeff_rot;

ConstantCoefficient coeff_vorticity;
FunctionCoefficient coeff_stream;
FunctionCoefficient coeff_stream_in;
ConstantCoefficient coeff_stream_closed_down;
ConstantCoefficient coeff_stream_closed_up;

//Interpolation objects
ParDiscreteLinearOperator gradient;
};

//Main class of the program
Expand Down Expand Up @@ -284,32 +277,26 @@ class Artic_sea{
ParMesh *pmesh;

FiniteElementCollection *fec_H1;
FiniteElementCollection *fec_ND;

ParFiniteElementSpace *fespace_H1;
ParFiniteElementSpace *fespace_ND;

Array<int> block_offsets_H1;

int dim;
double h_min;
int serial_refinements;
HYPRE_Int size_H1;
HYPRE_Int size_ND;

//System objects
ParGridFunction *temperature;
ParGridFunction *salinity;
ParGridFunction *relative_temperature;
ParGridFunction *phase;
ParGridFunction *vorticity;
ParGridFunction *stream;
ParGridFunction *velocity;
ParGridFunction *rvelocity;

BlockVector X;
BlockVector Y;
HypreParVector *Velocity;
HypreParVector *rVelocity;

//Solvers
Transport_Operator *transport_oper;
Expand Down Expand Up @@ -353,16 +340,16 @@ extern double InflowSalinity; //Salinity of the inflow
extern double NucleationSalinity; //Salinity of the nucleation point
extern double ZeroSalinity;

//Usefull position functions
//Usefull position functions~/Documents/brinicle/Brinicle
extern double r_f(const Vector &x); //Function for r
extern double r_inv_f(const Vector &x); //Function for 1/r
extern void zero_f(const Vector &x, Vector &f); //Function for 0 (vector)
extern void r_inv_hat_f(const Vector &x, Vector &f); //Function for (1/r)*r^ (r^ unitary vector)
extern void rot_f(const Vector &x, DenseMatrix &f); //Function for ( 0 1 )
// (-1 0 )

//Physical properties (in T,S)
extern double FusionPoint(const double T, const double S); //Fusion temperature at a given salinity
extern double FusionPoint(const double S); //Fusion temperature at a given salinity
extern double RelativeTemperature(const double T, const double S);
extern double Phase(const double T, const double S); //Phase indicator (1 for liquid and 0 for solid)
extern double HeatDiffusivity(const double T, const double S); //Coefficient for the diffusion term in the temperature equation
extern double SaltDiffusivity(const double T, const double S); //Coefficient for the diffusion term in the salinity equation
Expand Down
15 changes: 4 additions & 11 deletions Brinicle/code/run.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,13 +15,11 @@ Artic_sea::Artic_sea(Config config):
last(false),
vis_steps(config.vis_steps_max), vis_print(0),
pmesh(NULL),
fec_H1(NULL), fec_ND(NULL),
fespace_H1(NULL), fespace_ND(NULL),
fec_H1(NULL), fespace_H1(NULL),
block_offsets_H1(3),
temperature(NULL), salinity(NULL), phase(NULL),
temperature(NULL), salinity(NULL),
relative_temperature(NULL), phase(NULL),
vorticity(NULL), stream(NULL),
velocity(NULL), rvelocity(NULL),
Velocity(NULL), rVelocity(NULL),
transport_oper(NULL), flow_oper(NULL),
ode_solver(NULL), arkode(NULL),
paraview_out(NULL)
Expand Down Expand Up @@ -73,18 +71,13 @@ Flow_Operator::~Flow_Operator(){
Artic_sea::~Artic_sea(){
delete pmesh;
delete fec_H1;
delete fec_ND;
delete fespace_H1;
delete fespace_ND;
delete temperature;
delete salinity;
delete relative_temperature;
delete phase;
delete vorticity;
delete stream;
delete velocity;
delete rvelocity;
delete Velocity;
delete rVelocity;
delete transport_oper;
delete flow_oper;
delete ode_solver;
Expand Down
Loading

0 comments on commit c28e119

Please sign in to comment.