There are a lot of Mathematica packages for vector/tensor calculus: Tensorial, Advanced Tensor Analysis, Ricci, TensoriaCalc, grt, xAct are just a few mentions.

The main thing I don't like about these packages is that the majority of them (if not all) are not up to date with the current Mathematica version. Some of them haven't seen an update in decades. The other major drawback of these packages is the cumbersome notation, declaration and the fact that almost all are to be used with General Relativity in mind and use coordinate notation. The results are not really coordinate-free.

With this in mind, I have developed a "Package" that can deal with symbolic vectors/tensors in a coordinate-free form and the package is human-readable (some packages implementations in this matter are simply indecipherable).

In the first part, I'll present the code with simple explanations and the second part some simple examples.

These functions are needed to clear the OenOwnValues, DownValues, UpValues, SubValues. Mathematica has no built-in way of doing it.

(* Prevent evaluation *)
SetAttributes[{ClearOwnValues, ClearDownValues, ClearUpValues, ClearSubValues}, HoldFirst]

(* Always return true for commodity later *)
ClearOwnValues[var_Symbol]  := (OwnValues@var  = {}; True)
ClearDownValues[var_Symbol] := (DownValues@var = {}; True)
ClearUpValues[var_Symbol]   := (UpValues@var   = {}; True)
ClearSubValues[var_Symbol]  := (SubValues@var  = {}; True)

(* Delete Values of "f" if they match the input *)
ClearDownValues[expr:f_Symbol[___]]     := (DownValues@f = DeleteCases[DownValues@f, _?(!FreeQ[First@#, HoldPattern@expr] &)]; True)
ClearUpValues[expr:f_Symbol[___]]       := (UpValues@f   = DeleteCases[UpValues@f,   _?(!FreeQ[First@#, HoldPattern@expr] &)]; True)
ClearSubValues[expr:f_Symbol[___][___]] := (SubValues@f  = DeleteCases[SubValues@f,  _?(!FreeQ[First@#, HoldPattern@expr] &)]; True)

The Define-family functions are the main core of this package, they are used to define all kinds of relationship and can be easily expanded. You can define a Symbol "var" or a function "fun[var, ___]" to be a certain "type". "var" cannot have OwnValues, otherwise it will be evaluated. You can only define symbols. the function "fun" should be primarily used as script-like functions as Subscript, SuperHat, etc.

(* Prevent evaluation *)
SetAttributes[Define$Internal, HoldAll]

(* 
    Internal version of Define, all others versions are a call of this function.
    The OwnValues and/or DownValues are cleared immediately, this step is to avoid rule-parsing.
    The possible types are: Real, Imaginary, Tensor and Constant. For Tensor type additional parameters are needed: rank and dimension.
*)
Define$Internal[(var_Symbol /; ClearOwnValues@var) | (var:Except[Hold, fun_][head_Symbol /; ClearOwnValues@head, ___] /; ClearDownValues@var),
    type:"Real" | "Imaginary" | "Tensor" | "Constant",
    rank_Integer:2, dim_Integer:3] := Module[{tag},

    (* All expressions are defined as UpValues, assign it to the corresponding tag *)
    (* UpValues cannot be deeply nested, hence the need to assign it to the "tag" *)
    tag = If[Head@var === Symbol, var, fun];

    Which[
       type === "Real", (* Typical properties needed for real quantities *)
         Evaluate@tag /: Element[var, Reals] = True;
         Evaluate@tag /: Re[v:var] := v;
         Evaluate@tag /: Im[v:var] := 0;
         Evaluate@tag /: Conjugate[v:var] := v;
         Evaluate@tag /: Abs[v:var] := RealAbs@v;
       ,        
       type === "Imaginary", (* Typical properties needed for Imaginary quantities *)
         Evaluate@tag /: Element[var, Reals] = False;
         Evaluate@tag /: Re[v:var] = 0;
         Evaluate@tag /: Im[v:var] := v/I;
         Evaluate@tag /: Conjugate[v:var] := -v;
         Evaluate@tag /: Abs[v:var] := RealAbs@Im@v;
       ,        
       type === "Tensor", (* For compativility with Mathematica current Tensor-functions *)
         Evaluate@tag /: ArrayQ@var = rank != 0;
         Evaluate@tag /: TensorQ@var = rank != 0;
         Evaluate@tag /: MatrixQ@var = rank == 2;
         Evaluate@tag /: VectorQ@var = rank == 1;
         Evaluate@tag /: ListQ@var = rank != 0;
         Evaluate@tag /: ScalarQ@var = rank == 0;
         Evaluate@tag /: TensorRank@var = rank;
         Evaluate@tag /: TensorDimensions@var = ConstantArray[dim, {rank}];
         Evaluate@tag /: Element[var, Arrays@TensorDimensions@var] = True;
         Evaluate@tag /: Element[var, Matrices@{dim, dim}] = rank == 2;
         Evaluate@tag /: Element[var, Vectors@dim] = rank == 1;
       ,        
       type === "Constant", (* A constant has zero "derivative" *)
         Evaluate@tag /: ConstantQ@var = True;
         Evaluate@tag /: grad@var = 0;
         Evaluate@tag /: div@var = 0;
         Evaluate@tag /: curl@var = 0;
         Evaluate@tag /: DotNabla[_, var] = 0;
         Evaluate@tag /: D[var, __] = 0;  
         Evaluate@tag /: Dp[var, __] = 0; 
         Evaluate@tag /: Dt[var, ___] = 0;
         Evaluate@tag /: Delta[var] = 0;
       ,        
       True, $Failed]    
]

(* Assign more than one variable *)
Define$Internal[vars__ /; Length@{vars} > 1, type:"Real" | "Imaginary" | "Tensor" | "Constant", rank_Integer:2, dim_Integer:3] :=(
    Define$Internal[#, type, rank, dim] & /@ Hold /@ Hold@vars // ReleaseHold;) (* Hacky-way of passing Hold down *)
Define$Internal[Hold@var_, type:"Real" | "Imaginary" | "Tensor" | "Constant", rank_Integer:2, dim_Integer:3] := Define$Internal[var, type, rank, dim]

(* Main Define functions *)
SetAttributes[{DefineReal, DefineImaginary, DefineTensor, DefineConstant}, HoldAll]
DefineReal[vars__] := Define$Internal[vars, "Real"]
DefineImaginary[vars__] := Define$Internal[vars, "Imaginary"]
DefineTensor[vars__, rank_Integer:2, dim_Integer:3] := Define$Internal[vars, "Tensor", rank, dim]
DefineConstant[vars__] := Define$Internal[vars, "Constant"]

(* Define multiple things at once *)
SetAttributes[{DefineRealTensor, DefineConstantTensor, DefineRealConstantTensor}, HoldAll]
DefineRealTensor[vars__, rank_Integer:2, dim_Integer:3] := (DefineReal@vars; DefineTensor[vars, rank, dim];)
DefineConstantTensor[vars__, rank_Integer:2, dim_Integer:3] := (DefineConstant@vars; DefineTensor[vars, rank, dim];)
DefineRealConstantTensor[vars__, rank_Integer:2, dim_Integer:3] := (DefineReal@vars; DefineConstant@vars; DefineTensor[vars, rank, dim];)

Now it is possible to define tensorial variables and make them behave as tensor with current Mathematica implementation. Some built-in functions needed to be redefined to work with symbolic tensors. An example of this necessity is:

(* Define two tensors a and b *)
DefineTensor[a, b, 2]
TensorRank[2*a - 3*b] (* Return 2 *)
TensorQ[a] (* Return True *)
TensorQ[2*a] (* Return False *)

Mathematica function TensorQ don't know that a scalar times a tensor is a tensor. The following code is for refifinition:

Unprotect[TensorQ, VectorQ, TensorRank, Dot, Cross, TensorProduct]

(* Numbers are always scalar/constant. These functions are not built-in. *)
ScalarQ[a_?NumericQ]   := True
ConstantQ[a_?NumericQ] := True

(* Complexes *)
   TensorQ[(Re|Im|Conjugate)[a_]] := TensorQ@a
   VectorQ[(Re|Im|Conjugate)[a_]] := VectorQ@a
   ScalarQ[(Re|Im|Conjugate)[a_]] := ScalarQ@a
 ConstantQ[(Re|Im|Conjugate)[a_]] := ConstantQ@a
TensorRank[(Re|Im|Conjugate)[a_]] := TensorRank@a

(* Plus *)
TensorQ[(a_?TensorQ) + (b_?TensorQ)] := TensorRank@a === TensorRank@b
VectorQ[(a_?VectorQ) + (b_?VectorQ)] := True
ScalarQ[(a_?ScalarQ) + (b_?ScalarQ)] := True
ConstantQ[(a_?ConstantQ) + (b_?ConstantQ)] := True

(* Times *)
TensorQ[(a__?ScalarQ) * (b_?TensorQ)] := True
VectorQ[(a__?ScalarQ) * (b_?VectorQ)] := True
ScalarQ[(a__?ScalarQ) * (b_?ScalarQ)] := True
ConstantQ[(a_?ConstantQ /; ScalarQ@a) * (b_?ConstantQ)] := True
(* Pass scalars out of Dot and Cross, as is done in TensorProduct *)
Dot[a___, Times[b_, s__?ScalarQ], c___] := Times[s, Dot[a, b, c]]
Cross[a_, Times[b_, s__?ScalarQ]] := Times[s, Cross[a, b]]
Cross[Times[a_, s__?ScalarQ], b_] := Times[s, Cross[a, b]]

(* Dot *)
TensorQ[(a_?TensorQ) . (b_?TensorQ)] /; TensorRank@a + TensorRank@b - 2 >= 1 := True
VectorQ[(a_?TensorQ) . (b_?TensorQ)] /; TensorRank@a + TensorRank@b - 2 == 1 := True
ScalarQ[(a_?VectorQ) . (b_?VectorQ)] := True
ConstantQ[(a_?ConstantQ /; TensorQ@a) . (b_?ConstantQ /; TensorQ@b)] := True
(* Automatically evaluate to zero, as TensorProduct *)
Dot[a___, 0, b___] := 0

(* Cross *)
TensorQ[(a_?VectorQ) \[Cross] (b_?VectorQ)] := True
VectorQ[(a_?VectorQ) \[Cross] (b_?VectorQ)] := True
ConstantQ[(a_?ConstantQ /; VectorQ@a) \[Cross] (b_?ConstantQ /; VectorQ@b)] := True
(* Perpendicular vectors automatically evalute to zero *)
Cross[a_?VectorQ, a_?VectorQ] := 0
(* Automatically evaluate to zero, as TensorProduct *)
Cross[a___, 0, b___] := 0
(* Return single argument as Dot, Times and TensorProduct *)
Cross[a_] := a

(* Tensor Product *)
TensorQ[(a_?TensorQ) \[TensorProduct] (b_?TensorQ)] := True
ConstantQ[(a_?ConstantQ /; TensorQ@a) \[TensorProduct] (b_?ConstantQ /; TensorQ@b)] := True

(* Power *)
ScalarQ@Power[a_?ScalarQ, b_?ScalarQ] := True
ScalarQ[1/a_?ScalarQ] := True
ConstantQ@Power[a_?ConstantQ /; ScalarQ@a, b_?ConstantQ /; ScalarQ@b] := True
ConstantQ[1/a_?ConstantQ /; ScalarQ@a] := True

(* grad *)
grad[_?ConstantQ] := 0
TensorQ@grad[a_?ScalarQ] := True
VectorQ@grad[a_?ScalarQ] := True
TensorQ@grad[a_?TensorQ] := True
TensorRank@grad[a_?ScalarQ] := 1
TensorRank@grad[a_?TensorQ] := TensorRank@a + 1

(* div *)
div[_?ConstantQ] := 0
TensorQ@div[a_?TensorQ /; TensorRank@a >= 2] := True
VectorQ@div[a_?TensorQ /; TensorRank@a == 2] := True
ScalarQ@div[a_?VectorQ] := True
TensorRank@div[a_?TensorQ] := TensorRank@a - 1

(* curl *)
curl[_?ConstantQ] := 0
TensorQ@curl[a_?VectorQ] := True
VectorQ@curl[a_?VectorQ] := True
TensorRank@curl[a_?VectorQ] := 1

(* DotNabla *)
DotNabla[_, _?ConstantQ] := 0

(* Dp *)
Dp[_?ConstantQ, args__] := 0
TensorQ@Dp[a_, args__] := TensorQ@a
VectorQ@Dp[a_, args__] := VectorQ@a
ScalarQ@Dp[a_, args__] := ScalarQ@a
TensorRank@Dp[a_, args__] := TensorRank@a

(* Delta *)
Delta[_?ConstantQ] := 0
TensorQ@Delta[a_] := TensorQ@a
VectorQ@Delta[a_] := VectorQ@a
ScalarQ@Delta[a_] := ScalarQ@a
TensorRank@Delta[a_] := TensorRank@a

(* List *)
Dp[a_List, args__] := Dp[#, args] & /@ a

(* Don't assume anything is a scalar/constant *)
ScalarQ[a_]   := False
ConstantQ[a_] := False

Protect[TensorQ, VectorQ, TensorRank, Dot, Cross, TensorProduct]

Where we have defined the Tensor-functions: grad, div, curl which are self-explanatory; Dp is the partial derivative, Delta gives the variation of a quantity, somewhat related to Dp, and DotNabla is (for the lack of better name) the convective derivative.

For better print, we'll define the following notation:

(* Hacky-way to create parenthesis *)
MakeBoxes[Parenthesis[a_], _] := MakeBoxes[a.1][[1, 1]]
MakeBoxes[grad[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "grad", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Del]", #1} &)]
MakeBoxes[div[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "div", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Del]\[CenterDot]", #1} &)]
MakeBoxes[curl[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "curl", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Del]\[Cross]", #1} &)]
MakeBoxes[DotNabla[a_, b_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a, MakeBoxes@Parenthesis@b}, "DotNabla", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"(", #1, "\[CenterDot]\[Del]", ")", #2} &)]
MakeBoxes[Delta[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "Delta", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Delta]", #1} &)]

And now for the most important part of the code, the function ExpandDerivative, which as the name suggest, expand the derivative-like functions:

(* Expand Derivatives/Vectors/Tensors on expr and apply custom rules *)
ExpandDerivative[expr_, rules_:{}] := expr //. Flatten@{
    (* Custom Rules *)
    rules,

    (* Linearity *)
    (op:grad|div|curl|Delta|Inactive[grad]|Inactive[div]|Inactive[curl]|Inactive[Delta]|Re|Im|Conjugate)[a_ + b__] :> op@a + op[+b],
    (op:Dp|Inactive[Dp]|Sum|Inactive[Sum])[a_ + b__, arg__] :> op[a, arg] + op[+b, arg],
    (op:Times|Dot|TensorProduct|Cross|DotNabla|Inactive[DotNabla])[a___, b_ + c__, d___] :> op[a, b, d] + op[a, +c, d],
    (op:grad|div|curl|Delta|Inactive[grad]|Inactive[div]|Inactive[curl]|Inactive[Delta]|Re|Im|Conjugate)[(op\[CapitalSigma]:Sum|Inactive[Sum])[a_, args__]] :> op\[CapitalSigma][op@a, args],

    (* Sum *)
    (op:Sum|Inactive[Sum])[s_*a_, v_Symbol] /; FreeQ[s, v] :> s*op[a, v],
    (op:Sum|Inactive[Sum])[s_, v_Symbol] /; FreeQ[s, v] :> s*op[1, v],

    (* Complexes *)
    Conjugate@(op:Times|Dot|Cross|TensorProduct)[a_, b__] :> op[Conjugate@a, Conjugate@op@b], (* Pass Conjugate to child *)
    (op:grad|div|curl|Delta|Inactive[grad]|Inactive[div]|Inactive[curl]|Inactive[Delta])[(opC:Re|Im|Conjugate)[a_]] :> opC@op@a, (* Pass Conjugate/Re/Im to parent *)
    Dp[(op:Re|Im|Conjugate)[a_], v_] :> op@Dp[a, v], (* Pass Conjugate/Re/Im to parent *)

    (* Triple Product *)
    Cross[a_?VectorQ, Cross[b_?VectorQ, c_?VectorQ]] :> b*a.c - c*a.b,
    Cross[Cross[a_?VectorQ, b_?VectorQ], c_?VectorQ] :> b*a.c - a*b.c,

    (* Quadruple Product *)
    Dot[Cross[a_?VectorQ, b_?VectorQ], Cross[c_?VectorQ, d_?VectorQ]] :> (a.c)*(b.d) - (a.d)*(b.c),

    (* Second Derivatives *)
    div@curl[_?VectorQ] :> 0,
    curl@grad[_?ScalarQ | _?VectorQ] :> 0,

    (* grad *)
    grad[(s_?ScalarQ) * (b_)] :> s*grad@b + b\[TensorProduct]grad@s,
    grad[(a_?VectorQ) . (b_?VectorQ)] :> a\[Cross]curl@b + b\[Cross]curl@a + DotNabla[a, b] + DotNabla[b, a], (* Use physics form *)
    grad[(s_?ScalarQ) ^ (n_?ConstantQ)] :> n*s^(n-1)*grad@s,
    grad[(n_?ConstantQ /; ScalarQ@n) ^ (s_?ScalarQ)] :> n^s*Log[n]*grad@s,

    (* div *)
    div[(s_?ScalarQ) * (b_?TensorQ)] :> s*div@b + b.grad@s,
    div[(a_?VectorQ) \[TensorProduct] (b_?VectorQ)] :> DotNabla[b, a] + a*div@b,
    div[(a_?VectorQ) \[Cross] (b_?VectorQ)] :> b.curl@a - a.curl@b,

    (* curl *)
    curl[(s_?ScalarQ) * (b_?VectorQ)] :> grad[s]\[Cross]b + s*curl@b,
    curl[(a_?VectorQ) \[Cross] (b_?VectorQ)] :> div[a\[TensorProduct]b - b\[TensorProduct]a],

    (* DotNabla *)
    DotNabla[(s_?ScalarQ) * (b_?VectorQ), c_?VectorQ] :> s*DotNabla[b, c],
    DotNabla[a_?VectorQ, (\[Beta]_?ScalarQ)*(c_?VectorQ)] :> c*a.grad@\[Beta] + \[Beta]*DotNabla[a, c],

    (* Dp *)
    Dp[(op:Times|Dot|Cross|TensorProduct)[a_, b__], v_Symbol] :> op[Dp[a, v], b] + op[a, Dp[op@b, v]],
    Dp[Power[a_?ScalarQ, b_?ScalarQ], v_Symbol] :> Power[a, b-1]*b*Dp[a, v] + Power[a,b]*Log[a]*Dp[b, v],

    (* Delta *)
    Delta[(op:Times|Dot|Cross|TensorProduct)[a_, b__]] :> op[Delta@a, b] + op[a, Delta@op@b],
    Delta@Power[a_?ScalarQ, b_?ScalarQ] :> Power[a, b-1]*b*Delta[a] + Power[a,b]*Log[a]*Delta[b]
}

Some examples. Calculating the divergent of Maxwell Stress Tensor in vaccum:

Calculate the Einstein-Laub force density for linear dielectrics:

Testing Poynting theorem in vaccum (no sources):

Where the first argument is the quantity being "tested".

Many other uses are possible and is fairly easy to extended some definitions.